University of Southern California Dissertations and Theses

Dynamic modeling and simulation of flapping-wing micro air vehicles

(USC Thesis Other)

Dynamic modeling and simulation of flapping-wing micro air vehicles

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content Dynamic Modeling and Simulation of Flapping-Wing Micro Air Vehicles by Longlong Chang A Dissertation Presented to the FACULTY OF THE GRADUDATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulllment of the Requirements for the Degree DOCTOR OF PHILOSOPHY in Aerospace Engineering December 2019 Copyright 2019 Longlong Chang Dedication To my parents, and my elder brother Peifen Mao & Qingbin Chang, and Guo Chang ii Acknowledgments I would like to thank my advisor, Professor N estor O. P erez-Arancibia, for his guidance, encouragement, patience and support over the past six years. He provided me with the maximum freedom to explore various research topics, which broadened my horizons on interdisciplinary studies. As a researcher focusing on modeling and simulations, I really learned from N estor about the true nature of Engineering and essence of Science. I would also like to thank Professor Julian A. Domaradzki, Mitul Luhar and Professor Aiichiro Nakno for serving on my dissertation committee and their advice on my research and dissertation. Thanks to Professor Eva Kanso for serving on my qualifying committee and her critical comments on my research proposal. I want to thank all my labmates for their help, friendship and instructive academic interactions. Great thanks to Xiufeng Yang and Ariel Calder on for their collabo- rations on experiments for the validation of my modeling and simulation research. Thank Joey Zaoyuan Ge for his great friendship and experimental insights on non- linear dynamics. Thank Dr. Ying Chen for sharing of his solid knowledge on control and experiments. I also learned from him on how to plan and make progress. Thanks to Xuan-Truc Nguyen for her great help on my English writing and encouragements. Thank Emma Singer and Ke Coco Xu for discussions of research and great daily help. Also, thanks to my friends Xingtian Tao, Xize Wang and Wenhao Zhang for the shared life experience and happiness in Los Angeles. Thank USC and the Viterbi School for the resources supporting my research in the past years. Finally, great thanks to my parents, my elder brother and other family members for their long-term support, patience and love. iii Table of Contents Acknowledgments iii List of Figures viii List of Tables xvi Abstract xvii 1 Introduction 1 1.1 Background and motivation . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Contributions of the dissertation . . . . . . . . . . . . . . . . . . . . . 4 1.3 Organization of the dissertation . . . . . . . . . . . . . . . . . . . . . 5 2 Dynamics of a Single Wing in Passive Pitching 6 2.1 Problem description and assumptions . . . . . . . . . . . . . . . . . . 6 2.2 Dynamic modeling and integrated simulation . . . . . . . . . . . . . . 8 2.2.1 Equations of motion . . . . . . . . . . . . . . . . . . . . . . . 8 2.2.2 Integrated dynamics simulation with CFD solver . . . . . . . . 12 2.2.3 Quasi-steady models of a wing in passive pitching . . . . . . . 15 2.2.4 Pitching angle limit by hinge shoulders and wing deformation 21 2.2.5 Amended model for the dynamics of a wing in water . . . . . 23 iv 2.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.3.1 Experimental validation . . . . . . . . . . . . . . . . . . . . . 24 2.3.2 Applications of the integrated simulation with CFD . . . . . . 27 2.3.3 Adaptability of the revised QS models . . . . . . . . . . . . . 29 2.3.4 Integrated dynamics simulation with QS models . . . . . . . . 31 2.3.5 Pitching angle limit by hinge shoulders . . . . . . . . . . . . . 33 2.3.6 Integrated simulation for a wing immersed under water . . . . 34 3 Full Dynamics of the Entire Flapping-Wing MAV 35 3.1 Problem description . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.2 Lagrange's equations alternative . . . . . . . . . . . . . . . . . . . . . 38 3.3 Equations of motion . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.4 Singularity-free and time-averaged formulations . . . . . . . . . . . . 47 3.4.1 Formulation with dual sets of Euler angles and quaternions . . 47 3.4.2 Time-averaged models with mass-less wing assumption . . . . 50 3.5 Integrated dynamics simulation with CFD solver and controller . . . 53 3.6 Quasi-steady models for a MAV in free ight . . . . . . . . . . . . . . 58 3.6.1 Wing with arbitrary root kinematics . . . . . . . . . . . . . . 58 3.6.2 Simplied cylindrical body . . . . . . . . . . . . . . . . . . . . 60 3.7 Preliminary experimental validation . . . . . . . . . . . . . . . . . . . 63 3.7.1 Options for experimental validation . . . . . . . . . . . . . . . 63 3.7.2 Experimental setups for underwater and in-air tests . . . . . . 65 3.7.3 Modeling modications for the underwater tests . . . . . . . . 65 3.8 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 3.8.1 Integrated simulation with CFD solver and controller . . . . . 68 v 3.8.2 Adaptability of the quasi-steady models for the wing and body in free ight . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 3.8.3 Full dynamics simulation with quasi-steady models . . . . . . 74 3.8.4 Experimental validation of the full dynamics model with un- derwater tests . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 3.8.5 Dynamic behavior analysis for MAV II ight in the air . . . . 82 3.8.6 Eects of the added mass correction, initial conditions and hinge stiness . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 3.8.7 Towards new designs of apping-wing MAVs . . . . . . . . . . 85 4 Unied CFD Solver Prototype for Moving Boundary Problems 88 4.1 Implementation strategies for the unied CFD solve prototype . . . . 90 4.1.1 Arbitrary Lagrangian-Eulerian approach and overset grid method 90 4.1.2 Overset grid assembling and discretization schemes for the seg- regated solver . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 4.1.3 Solver implementation and parallelization strategies . . . . . . 97 4.2 Fast overset grid assembling and automatic parallel multigrid generation 99 4.2.1 Search-based fast overset grid assembling algorithm . . . . . . 99 4.2.2 Hybrid parallel automatic multigrid generation . . . . . . . . . 104 4.3 Smoothness and consistency enhancements . . . . . . . . . . . . . . . 109 4.3.1 High-order face interpolation for non-regular mesh . . . . . . . 109 4.3.2 Smoothness and consistency of FVM for small t . . . . . . . 111 4.3.3 Coupling strategies for interpolation cell of the overset grids . 114 4.4 Testing examples for OgFSI . . . . . . . . . . . . . . . . . . . . . . . 115 4.4.1 Overset grid assemblies with search-based algorithm . . . . . . 116 4.4.2 Cell index re-numbering and multigrid generation . . . . . . . 116 vi 4.4.3 Finite volume solver based on overset grid . . . . . . . . . . . 118 4.4.4 Results for moving cylinder problems . . . . . . . . . . . . . . 121 4.4.5 Results for moving wings . . . . . . . . . . . . . . . . . . . . . 121 5 Conclusions and Future Work 123 Bibliography 126 Appendices 135 Appendix A. Matrices and derivatives . . . . . . . . . . . . . . . . . . . . . 135 Appendix B. Parameters of the wing and body . . . . . . . . . . . . . . . . 139 Appendix C. File locking for parallel implementation of User Code . . . . . 141 vii List of Figures 1.1 Study of apping wing yers at insect-sized scale. (a) Theoretical analysis and high lift generation mechanisms [1{4]. (b) Quantitative experimental measurements of aerodynamic forces [5{8]. (c) Flight dynamics [9{12]. (d) Robotic design, fabrication and control [13{19]. 2 2.1 Conceptual design of the studied apping-wing MAVs and supporting information for dynamic modeling. (a) Dipteran inspired MAV with a non-uniform cylindrical body and dual wings. (b) A video frame of the wing status for a static test in the air. (c) Computing process under the weakly-coupled FSI assumption. . . . . . . . . . . . . . . . . . . . 6 2.2 Descriptions of wing kinematics and exure hinge geometry. (a) Wing rotations with respect to a xed point and frames of reference. (b) Simple beam model for the exure hinge. . . . . . . . . . . . . . . . . 9 2.3 Integrated dynamics simulation and overset grids. (a) Data ow for the integrated FSI simulation. (b) Illustration of overset grid assem- bly with elliptic wing airfoil. (c) Assembly of 3D mesh blocks. (d) Topology for body-tted mesh of the wing with a regular interface. (e) Topology for body-tted mesh of the wing with a non-regular interface. 13 viii 2.4 Quasi-steady components of aerodynamic forces on one wing section for hovering ight. (a) Translational component due to angle of attack and incident velocity. (b) The rst rotational component to satisfy Kutta condition. (c) The second rotational component due to damping. (d) Added mass component due to wing acceleration. . . . . . . . . . . . 15 2.5 Geometric limit by hinge shoulders and correction methods. (a) Geo- metric limit applied by hinge shoulders. (b) Geometric limiter method. (c) Enlarged stiness and damping method. . . . . . . . . . . . . . . 21 2.6 Experimental setup and measured kinematics from [7] for validation of the single wing dynamics. (a) Illustration of the experimental setup. (b) Planar shape of the wing in experiment and simulations. (c) Ex- perimentally measured apping angle, pitching angle and deviation angle in [7]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.7 Model validation for the dynamics of passive wing pitching at 100 Hz. (a) Instantaneous pitching angles comparison. (b) Vertical in- ertial force computed from simulated and measured . (c)-(d) Total vertical force with dierent ltering methods for measured . FS rep- resents Fourier series tting and ZP denotes zero-phase ltering. . . . 25 2.8 Parametric studies of hinge stiness and wing planar shape. (a)-(b) Pitching angles and vertical aerodynamic forces for various hinge sti- ness. (c)-(d) Pitching angles and vertical aerodynamic forces for a Eristalis-tenax-like wing and a honeybee-like wing for a hinge stiness k h = 3. The unit for hinge stiness is Nm. . . . . . . . . . . . . . . 28 2.9 CFD validation of quasi-steady models for a single wing in hovering. (a)-(c) Instantaneous lift, drag and pitching torque for t = 0:5. (d)-(f) Instantaneous lift, drag and pitching torque for t =. . . . 30 ix 2.10 Quasi-steady model calibration for a single wing in passive pitching. (a)-(d) are instantaneous pitching angle, pitching moment, lift and drag, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.11 Pitching angle limit with the geometric limiter method. . . . . . . . . 33 2.12 Converged integrated simulations for underwater tests at 13 Hz. (a) Wing planar shape. (b) Converged pitching angle proles. . . . . . . 34 3.1 Free ight and kinematic description of the MAV. (a) An illustration of the free ight. (b) Frames of reference. (c) Flapping and pitching angles with respect to the body. . . . . . . . . . . . . . . . . . . . . . 36 3.2 The dual sets of Euler angles. (a) Coordinate frames describing the dual sets of Euler angles. (b) The relative rotation for the calculation of the roll angle 0 b . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.3 Computation work ows for various MAV dynamics models. . . . . . 50 3.4 Integrated dynamics simulations with CFD solver and controller. (a) Integrated simulation for the entire MAV. (b) Mesh topology for the CFD simulation of the entire MAV. (c) One slice for the assembled grids. (d) Simulation with controller included. . . . . . . . . . . . . . 54 3.5 Parallel implementation of the integrated simulation with STAR-CCM+ and MATLAB. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 3.6 Blade element analyses for the wing and the body in free ight. (a) A wing with an arbitrary wing root kinematics. (b) Body motion and wing induced air ow. . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 3.7 The 70-mg insect-sized robotic MAV and related experimental settings. (a) CAD model for the MAV II. (b) Front view of the open-loop test in the air. (c) Front view of the open-loop test in water tank. . . . . . 63 x 3.8 Integrated dynamics simulations with CFD solver for MAV I. (a) In- stantaneous apping angles. (b) Planar position and attitude. (c) Body pitch angle. Note: the open-loop simulation terminated around 0.19 s because the moving mesh reaches the boundary of the back- ground mesh. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 3.9 Flow elds for the open-loop ight of the MAV. (a)-(b) Pressure and velocity magnitude elds for the longitudinal plane of the body. (c)- (d) Pressure and velocity magnitude elds for the wing section with maximum chord length. . . . . . . . . . . . . . . . . . . . . . . . . . 70 3.10 CFD validation of the quasi-steady model for the case with incoming air ow. (a)-(c) Lift, drag, pitching moment for U = 0:5m/s. (d)-(f) Lift, drag, pitching moment for U = 1:0m/s. . . . . . . . . . . . . . . 72 3.11 CFD validation of the quasi-steady model for the body force resulting from the wing induced air ow. . . . . . . . . . . . . . . . . . . . . . . 73 3.12 Comparisons of the body states between the CFD and QS results for MAV I. (a) Planar coordinates of the body centroid for 0 = 0. (b) Body pitch angle for 0 = 0. (c) Body pitch angle for 0 = 0:1. . . . . 75 3.13 Comparisons of the instantaneous forces on the left wing between the CFD and QS results. (a) Instantaneous vertical force in the inertial frame. (b) Instantaneous horizontal force in the inertial frame. (c) Pitching moment of the wing about the hinge axis in wing-xed frame. 76 3.14 Full dynamics model validation with the underwater test for b 0 = 0:05. (a) Dynamic behavior comparison between the experimental (left) and simulated (right) results. See the clip C4 in supplementary video S1.mp4 for more frames. (b) Instantaneous planar coordinates of the MAV body. (c) Instantaneous pitch angle of the body. . . . . . 79 xi 3.15 Full dynamics model validation with the underwater test for b 0 = 0:0. (a) Dynamic behavior comparison between the experimental and simu- lated results. See the clip C5 in supplementary video S1.mp4 for more frames. (b) Instantaneous planar coordinates of the MAV body. (c) Instantaneous pitch angle of the body. . . . . . . . . . . . . . . . . . 80 3.16 Comparisons of the experimental ying behaviors with the simulated results. (a) The baseline jump test with nominal bias b 0 = 0:0. (b) The test with large bias correction. (c) The test a small bias correction. 83 3.17 Eects of added mass correction, initial conditions and hinge sti- ness. (a) Sensitivity analysis on the added correction coecient. (b) Dynamic response with/without the ground interaction model for dif- ferent initial conditions. (c) Dynamic responses for dierent values of hinge stiness without the collision model and ground interaction model. 85 3.18 New design of apping-wing MAVs inspired or potentially investigated by the studies of full dynamics models. (a) Photo sequence for the ight of Bee + , adapted from [18]. (b) Pitch oscillation of the body for low-frequency butter y-like ying. . . . . . . . . . . . . . . . . . . . . 86 4.1 Common Lagrangian approaches for mesh boundary updating. (a) Immersed boundary method. (b) Dynamic mesh method. (c) Sliding grid method. (d) Overset grid method. (e) Boundary approximation with region domain mesh. . . . . . . . . . . . . . . . . . . . . . . . . 91 4.2 Overset grid assembly. (a) Original mesh blocks. (b) Articial bound- ary after hole-cutting. (c) Interpolation cell and stencil. . . . . . . . . 93 xii 4.3 Finite volume discretization for a mesh cell. (a) Mesh neighboring information. (b) Gradient estimation with deferred correction, adapted from [95]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 4.4 Technical strategies for the implementation of OgFSI. . . . . . . . . . 97 4.5 Parallel implementation with global indexing. . . . . . . . . . . . . . 98 4.6 Distance-based hole cutting, neighbor-based linear host searching adapted from [82] and index re-numbering. . . . . . . . . . . . . . . . . . . . 99 4.7 Mesh topologies for overset grid assembling. (a) One layer with mul- tiple moving objects. (b) Two layers. (c) Partially overlapping mesh blocks. (d) On mesh block partially outside background mesh. . . . . 101 4.8 Search-based fast hole-cutting algorithm.(a)-(c) Hole cutting for a bound- ary with continuous curvature. (d)-(f) Hole cutting for a boundary with non-continuous curvature. . . . . . . . . . . . . . . . . . . . . . 102 4.9 Matrix topology and AMG process for the overset grid-based solver. . 104 4.10 Automatic agglomeration algorithm for serial and parallel generators. 106 4.11 Non-smoothness of velocity eld for non-regular mesh. (a) Non-regular mesh for channel ow. (b) Result of velocity eld from OgFSI-old. (c) Result from STAR-CCM+ 12.02. (d) Result from Fluent 15.2. (e) Result from OgFSI with a high-order interpolation scheme. . . . . . . 110 4.12 High-order interpolation for the updating of face variables considering the non-regular shape of the cell. . . . . . . . . . . . . . . . . . . . . 111 4.13 Non-smoothness of the pressure eld at the 2nd time step. (a) Multi- block mesh for the cylinder ow. (b) OgFSI-old at t=0.1 s. (c) OgFSI-old at t=0.1 ms. (d) STAR-CCM+ at t=0.1 ms. (e) Fluent at t=0.1 ms. (f) OgFSI at t=0.1 ms. . . . . . . . . . . . . . . . . 112 xiii 4.14 Non-smoothness of the velocity and pressure elds for overset grid solu- tions at small time step t = 0:1 ms. (a) Velocity eld by Fluent. (b) Pressure eld by STAR-CCM+. (c) Smooth velocity eld by OgFSI. (d) Smooth pressure eld by OgFSI. . . . . . . . . . . . . . . . . . . 114 4.15 Strong coupling vs loose coupling for solving pressure correction. (a) is re-plotted from Fig. 4.2(c). . . . . . . . . . . . . . . . . . . . . . . 115 4.16 Overset grid assembly with the search-based fast algorithm. . . . . . 116 4.17 Cell index re-numbering and coarsening for single block mesh. . . . . 117 4.18 Examples of the automatic multigrid generator for overset grids. (a)- (b) Automatic agglomeration for serial and parallel coarsening. (c)-(d) Multigrid generation for multiple mesh blocks and 3D meshes. . . . . 118 4.19 x velocity eld for the lid-driven cavity problem. (a) Single structured quadrilateral mesh block with 100 100 cells. (b) Single unstructured triangular mesh block. (c) x velocity eld and mesh assembling for OgFSI and STAR-CCM+ with overset grids. (d) FVM-based IBM method using the overset grid hole-cutting. . . . . . . . . . . . . . . . 119 4.20 Results for the moving rigid overset grid solvers for unsteady simulation in 150 s. (a) Rigid cylinder with STAR-CCM+. (b) Rigid cylinder with OgFSI. (c) Deformable cylinder with OgFSI. (d) Mesh deformation for the cylinder with OgFSI. . . . . . . . . . . . . . . . . . . . . . . . . . 120 4.21 x velocity elds for the moving wings with OgFSI. (a) Rigid 2D elliptic wing with prescribed horizontal velocity and pitching angular velocity. (b) Rigid 3D rectangular wing. (c)-(d) Vertical elliptic wing in hori- zontal motion with enlarged minor axis. (e)-(f) Moving elliptic wing with local bending. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 xiv A.1 Geometry and diagram of the wings and body for the two MAVs. (a)- (b) MAV I. (c)-(d) MAV II. . . . . . . . . . . . . . . . . . . . . . . . 139 xv List of Tables 2.1 Parameters of the wing for the benchmark passive simulation. . . . . 26 2.2 Cycle-averaged lift for various hinge stiness. . . . . . . . . . . . . . 28 2.3 Parameter sets for the quasi-steady models of a single wing. . . . . . 29 2.4 Parameter sets for the quasi-steady model of a wing in passive pitching. 32 3.1 Parameters for the quasi-steady model of the body. . . . . . . . . . . 74 3.2 Parameter set for the QS model of the wing in underwater test. . . . 78 3.3 Parameter set for the QS model of the body in underwater test. . . . 78 3.4 Parameter set for the model corrections in simulating underwater tests. 81 4.1 Cell index re-numbering algorithm for each mesh block. . . . . . . . . 100 4.2 Proposed fast search-based hole-cutting algorithm. . . . . . . . . . . . 103 4.3 Eciency testing for the hole-cutting algorithms. . . . . . . . . . . . 104 4.4 Serial search algorithm for the coarser layer of grid. . . . . . . . . . . 107 4.5 Parallel search algorithm for the coarser layer of grid. . . . . . . . . . 108 A.1 Geometric parameters for the wings. . . . . . . . . . . . . . . . . . . 139 A.2 Inertial parameters for the wings. . . . . . . . . . . . . . . . . . . . . 139 A.3 Geometric parameters for the bodies. . . . . . . . . . . . . . . . . . . 140 A.4 Inertial parameters for the bodies. . . . . . . . . . . . . . . . . . . . . 140 xvi Abstract Birds and insects inspired us to develop ying machines, such as airplanes and heli- copters. However, the development of articial apping-wing yers falls far behind that of other aircraft, which is possibly due to the complexities of the associated unsteady aerodynamics, bidirectional actuation, stability and control, transmission design, etc. To aid the practical design of insect-sized robotic micro air vehicles (MAVs), we aimed to develop a set of tools to analyze the unsteady aerodynamics and ight dynamics of apping-wing yers. The robotic design of MAVs with wings connected to the body through exure hinges is able to reduce the number of actuation degrees of freedom and the complexity of transmission design. In this design, each wing is partially actuated and the pitching motion of the wing is dynamically determined under the interactions among inertial, hinge elastic, and aerodynamic forces. Both the wing motion itself and the dynamics of the entire MAV dene uid-structure interaction (FSI) procedures. For estimating the aerodynamic forces on a single wing in passive pitching, the dynamics model of one wing is rst derived and simulated by integrating the equations of motion with a proprietary computational uid dynamics (CFD) solver. To achieve physical insights and quick force estimation, revised empirical formulas are proposed to extend the adaptability and improve the accuracy of the existing quasi-steady (QS) models. Preliminary models for the pitching angle limit applied by hinge geometry are proposed and simulated. The validity of the weakly decoupled FSI assumption is discussed with the simulation of a single wing immersed in water. The full dynamics model of the entire MAV is formulated using our proposed alternative form of the Lagrange's equations. In this alternative form, repeated in- xvii ertial terms are removed to enable the manual derivation of compact equations of motion for the MAV to achieve ecient numerical computations. The integrated dynamics simulation with the proprietary CFD solver provides us with high-delity numerical predictions for the dynamic behavior of the apping-wing MAVs in open- loop or closed-loop simulations. The computationally ecient full dynamics model is formulated with QS formulas. Simulation results are compared with the results by the integrated simulation with CFD solver. The full dynamics model with QS formulas is experimentally validated using both the underwater and in-air jump tests for quantitative and qualitative purposes, respectively. Highly satisfactory matching between the simulated and experimental results is achieved. The full dynamics model is proved to be eective in the dynamic behavior prediction for the robotic MAVs. To compensate for the incapability of the proprietary CFD solver in exibility and computing eciency, a unied CFD solver based on the overset grid method is proposed and implemented for moving boundary problems at low Reynolds numbers. To achieve ecient and accurate simulations with the segregated overset grid solver, a fast search-based hole-cutting algorithm, a hybrid automatic multigrid generation process, and several consistency-improved discretization schemes are proposed. Fi- nally, testing examples are provided to demonstrate the capability of the in-house solver on unsteady and moving boundary problems. xviii Chapter 1 Introduction 1.1 Background and motivation Human beings have realized the dream of ying in the sky with the splendid in- ventions of airplanes, helicopters, airships and rockets. In the process of developing articial ying machines, we were inspired by the wings of birds and insects, and developed the theory of wings for aircraft design. However, engines or motors with unidirectional rotating components or air injection are the dominating means of actu- ation generation for articial yers. Practical diculties, including highly nonlinear unsteady aerodynamics, mechanical complexities, power eciency and materials, pre- vented us from developing mature prototypes of apping-wing air vehicles. To achieve a better understanding of apping-wing ight, during the past decades, the study of apping-wing yers at insect size has been widely conducted in several categories as illustrated in Fig. 1 (a)-(d). Signicant progress has been made in (1) theoretical analysis and high lift mechanisms of unsteady aerodynamic force generation on a wing at low Reynolds number, (2) quantitative experimental measurements of aerodynamic forces on wings or yers, (3) ight dynamics investigation of the entire yers, and (4) robotic design, fabrication and control. Among the robotic MAVs in Fig. 1(d), the development of real insect-sized proto- types [17,18] was achieved by employing a novel fabrication methodology: the smart 1 Taylor and Thomas (2003) Sun and Xiong (2005) Chang and Wang (2014) Orlowski and Girard (2011) Weis-Fogh and Jensen (1956) Weis-Fogh (1973) Ellington, et al(1996) Sane (2003) (b) (c) (a) Weis-Fogh (1956) Dickinson, et al. (1999) Whitney and Wood (2010) Singer, et al. (AMSL, 2019) (d) M. Karasek, et al. (2018) X. Yang, et al. (USC, 2019) Z. Tu, et al. (2019) S. B. Fuller (2019) Roshanbin, et al. (2019) S. Deng, et al. (2014) Festo (2013) R. J. Wood (Harvard, 2008) Fig. 1.1. Study of apping wing yers at insect-sized scale. (a) Theoretical analysis and high lift generation mechanisms [1{4]. (b) Quantitative experimental measure- ments of aerodynamic forces [5{8]. (c) Flight dynamics [9{12]. (d) Robotic design, fabrication and control [13{19]. composite microstructures (SCMs). The SCM technology enables a compact design: each wing aps actively under the actuation of transmission, and pitches passively around the spanwise axis of a exure hinge. Thus, the wings are not fully actively actuated, and the pitching motion of the wing about the exure hinge is referred to as the passive rotation mechanism in this study. For each wing, since only one degree of freedom (DOF) is directly driven by the transmission, the passive rotation mechanism can eectively reduce the complexity of the transmission design and the total weight of the entire MAV. It is, therefore, a crucial step toward the successful development of MAVs at the insect scale. 2 The ight of a MAV in space is a FSI process. The eciency and eectiveness of a MAV heavily depend on the performance of aerodynamic force generation on wings. For a wing with exure hinge, the pitching motion has its own dynamics, which is an interactive process between the inertial and aerodynamic forces. Thus, the passive rotation mechanism brings more challenges to the experimental and theoretical studies of apping wings. For example, the scaled-up experiments based on Reynolds number (Re) matching for aerodynamic force measurement in [6, 20, 21] cannot be directly performed for the passive rotation case. Besides the scaling of geometric dimensions, hinge stiness and wing inertial properties must also be correctly scaled. After wing inertia is scaled-up, the inertial force of a wing is more dicult to be separated from the measured total force. For numerical studies, even the aerodynamic forces estimation of a single wing is a numerical FSI problem [7, 22, 23], which is more challenging to obtain a closed-form solution or approximation than that of a wing with fully prescribed wing kinematics. Theoretical studies of apping wings with passive rotation mechanisms have been primarily carried out for the aerodynamic analysis of a single wing with its rotating point xed on the ground or hovering in the air [7,22{25]. To answer questions such as guidelines for installing wings onto the airframe or the quantitative eects of wing inertia on the open-loop dynamic response of a MAV, the ight dynamics model of the entire air vehicle must be investigated. Most research on the dynamics of the apping- wing ight is for the yers with fully prescribed wing kinematics or a cycle-averaged pitching angle as exemplied in Fig. 1(c) [9, 11, 12, 26{33]. Only limited works were directly or indirectly carried out for the yer with passive hinges using the Newtonian methods [12,34]. Therefore, it is interesting and useful to investigate the dynamics of this kind of MAV in a more generalized multibody dynamics approach. To achieve a high-delity estimation of aerodynamic forces acting on the MAV, the computational uid dynamics solver with moving boundary capability was employed [23,25,27,35,36]. However, the availability of the CFD solver can be a barrier for other researchers 3 to work on related topics. A rapid and unied CFD solver for moving boundary problems at low Reynolds numbers is necessary to complete the chain of analysis tools for apping-wing studies. During my research in the past Ph.D. years, I aimed to formulate dynamic models and develop a set of tools for analyzing the dynamic behaviors of the apping-wing MAVs to aid our mechanical and controller designs of the real robotic systems. 1.2 Contributions of the dissertation The main contributions of this dissertation are summarized as follows: Dynamics model for a single wing in passive pitching was derived and simulated by incorporating the proprietary CFD solver and revised quasi-steady models into the equation of pitching motion. On the revisions of the existing quasi-steady models, enhanced adaptability and accuracy were achieved and validated with CFD simula- tions. The preliminary models for the pitching angle limit by hinge geometry were proposed and applied in the full dynamics simulation. The weakly coupled equation of motion for a single wing in the air was modied and extended to the simulation of a single wing immersed in water. An alternative form of Lagrange's equations was derived through removing re- peated terms, which provides us with a numerically ecient formulation for the com- putation of inertial forces. With this alternative form, the full dynamics model of the entire MAV was formulated and integrated with the proprietary CFD solver and QS models. The underwater tests provided a preliminary quantitative validation for the proposed full dynamics model. Simulations fo jump tests in the air and other design applications showed the capability of the full dynamics model on the prediction of the dynamic behavior of apping-wing MAVs. A unied CFD solver prototype was proposed and implemented for moving boundary problems. A fast overset grid hole-cutting algorithm, a hybrid parallel 4 multigrid generation process and several consistency enhanced discretization schemes for the segregated solver were proposed and demonstrated with numerical tests. 1.3 Organization of the dissertation The dissertation is organized as follows: The dynamics of a single wing in passive pitching motion is discussed in Chapter 2. The fundamental assumptions and sim- ulation approaches for this dissertation are presented in detail. In Chapter 3, the full dynamics model of the entire MAV is derived and numerically simulated with CFD solver and quasi-steady models. Preliminary experimental validations in water and air are provided to demonstrate the capability of the full dynamics model. The singularity-free and time-averaged models are also discussed. Chapter 4 describes the implementation and enhanced algorithms of the unied in-house CFD solver based on the moving overset grid method. Finally, conclusions and future work are discussed in Chapter 5. 5 Chapter 2 Dynamics of a Single Wing in Passive Pitch- ing 2.1 Problem description and assumptions The insect-sized apping-wing MAVs studied in our research, as depicted in Fig. 2.1(a), have the basic robotic design developed in [17,18,37]. In this robotic design, wings are fabricated with polyester lm and reinforced with carbon ber frames. The airframe and transmission are constructed with carbon bers based on the SCMs methodol- ogy [17]. From the conceptual design perspective, the airframe is simplied as in a Body Wing Flexure hinge Polyester film End effector Carbon fiber frame (b) (c) (a) Moving Boundaries Fluid-Structure Interaction Aerodynamic Forces Equation of Motion Flow Equations Fig. 2.1. Conceptual design of the studied apping-wing MAVs and supporting information for dynamic modeling. (a) Dipteran inspired MAV with a non-uniform cylindrical body and dual wings. (b) A video frame of the wing status for a static test in the air. (c) Computing process under the weakly-coupled FSI assumption. 6 cylindrical rigid shape. The use of this simplication will be explained in the next chapter. As discussed in Sec. 1.1, one fundamental feature of this kind of MAVs is that the wings are connected to the end eectors of the transmission through exure hinges, that is, the passive rotation mechanism. During the ight of the MAV, the instantaneous pitching angle of each wing is dynamically determined by the inter- actions between wing inertia, hinge elasticity, and aerodynamic forces on the wing. Compared to a MAV with fully actively actuated wings, the passive rotation mech- anism brings more challenges to the precise description of the MAV's states and the prediction of the MAV's dynamic behavior. To approach the dynamic modeling of the considered MAVs, several important assumptions are employed in our discussion. First, the real inputs to the MAV are the electrical signals sent to the piezoelectric actuators. Wings ap under the driv- ing torque by the end eectors of the actuator-transmission system. The actuation system preserves its own dynamic process. In [38], a time-invariant linearized model was proved to be accurate enough for the mapping from the input electrical sig- nal to the apping angle of the wing. The inclusion of the actuator-transmission system brings more parameters and uncertainties to the wing dynamics. To avoid the modeling of the actuation system, the apping angles are directly described in prescribed/controller-output proles as the input signal to the MAV instead of the electrical signal to actuators. The prescribed apping angles are assumed to be al- ways achievable for the considered MAV. Second, the wings are reinforced with carbon bers to enhance the stiness of the structure. From the video frame of our static test shown in Fig. 1.1(b) and related observations in [7], for properly selected wing kinematics (e.g., apping frequency around 100 Hz), this kind of wing can maintain its planar shape without signicant deformation except for some upstroke-downstroke transition frames. In this study, wings are treated as rigid thin plates to reduce the complexity of dynamics modeling. Third, as the functional component for passive pitching, the exure hinge is fabricated with Kapton lm and modeled as torsional 7 springs with constant stiness in Sec. 2.2. With the above three assumptions, we are able to derive the equations of motion for the wing as a discrete and lumped system. However, the dynamics of the wing is a FSI process, and the equations of motion and Navier-Stokes equations cannot be solved simultaneously in a simple method. To make use of existing analysis tools, the wing dynamics is assumed as a weakly coupled FSI process. We can solve the equations of motion for the wing and the governing equations for the ow separately and sequentially. The equation of motion of the wing provides a moving boundary for the ow eld and the ow eld equations generate aerodynamic forces on the wing. These assumptions are also applied to the dynamic modeling of the entire MAV in the next chapter. With the above assumptions, we can minimize the number of parameters to de- scribe the system and reduce the uncertainties from wing exibility and actuation system. Also, with the weakly-coupled FSI assumption, we can access to more avail- able tools easily. 2.2 Dynamic modeling and integrated simulation 2.2.1 Equations of motion In this chapter, the wing is assumed in ground test or hovering status, that is, the root of the wing is xed in space as depicted in Fig. 2.2 (a). Since the wing is treated as a rigid plate and the exure hinge is modeled as a torsional spring, the Lagrange's equations can be employed to derive the equations of motion for this kind of discrete system. For each generalized coordinate, Lagrange's equations provide us with scalar forms for the generalized inertial and potential forces, but the original form of Lagrange's equations lacks an explicit formula for the non-conservative forces [39]. 8 (a) (b) Fig. 2.2. Descriptions of wing kinematics and exure hinge geometry. (a) Wing rotations with respect to a xed point and frames of reference. (b) Simple beam model for the exure hinge. The Lagrange's equations in Kane's formulation are employed [40], d dt @L @ _ q j @L @q j = Nc X i=1 f i @v i @ _ q j + i @! i @ _ q j ;j = 1; 2;:::;N q : (2.1) Note that, on the left-hand side, L = TV is the Lagrangian for the rigid body system, T is the kinetic energy, V is the potential energy, q j is one generalized co- ordinate, and N q is the number of the generalized coordinates. The coordinates q j s are independent of each other. On the right-hand side, N c is the number of compo- nents in the system, v i is the translational velocity of the center of mass of the ith component, ! i is the angular velocity of the ith component. f i and i denote the non-conservative force and moment acting on the center of mass of theith component, respectively [40]. To describe the motion of the wing, two frames of reference are dened as shown in Fig. 2.2(a). F 0 fo 0 ;x 0 ;y 0 ;z 0 g denotes a right-handed inertial frame, which has its origin coinciding with the wing rotation point, itsx 0 axis vertically upwards and its y 0 axis pointing outside along the initial position of the hinge axis. F w fo w ;x;y;zg denotes a right-handed wing-xed frame with the same origin as that of F 0 (o 0 is displaced for gure clarity), with its axisy extending along the instantaneous position of the hinge axis and itsz axis normal to the wing plane. The attitude of the wing is described by three consecutive Euler angles betweenF w andF 0 , that is, apping angle 9 aboutx 0 axis, pitching angle about the carriedy axis and stroke plane deviation angle about thez axis. For the design in Fig. 2.2(a), there is no specic mechanism for the deviation motion except for some small elastic eects. For simplicity, we neglect the deviation motion ( 0 and _ 0) and only consider the apping and pithing angles for the dynamic modeling of the wing. As depicted in Fig. 2.2(b), the exure hinge is treated as a simple single-axis bending beam as the model in [7]. Since the exure hinge is the most crucial part of the passive rotation mechanism, we want to explain the modeling of the hinge in detail. The hinge has one face xed to the driving spar and deforms according to a circular arc. Thus, from Hooke's law and the denition of strain, the curvature of the neutral surface can be computed as = h (EI h ) 1 , where h is the externally applied moment to the beam, E is the elastic modulus of the hinge material and I h is the second moment of area of the hinge transverse section (I h = 12 1 w h t 3 h with width w h and thicknesst h ). For the deformation in Fig. 2.2(b) (circular arc), the de ection angle of the exure hinge is =l h , where l h is the length of the exure hinge. The equivalent torsional stiness of the exure hinge is computed as k h = EI h l h = Ew h t 3 h 12l h : (2.2) To determine the exact instantaneous attitude of the wing, the position of the equiv- alent rotation axis of the wing must be estimated. From Fig. 2.2(b), the critical parameter that denes the position of the wing's rotation axis, the distance h o , is h o 1 sin 2 = l h sin 2 : (2.3) Taking into account the maximum possible geometric de ection of the hinge and the relatively small value ofl h compared to the wing chord, the rotation axis shift during wing pitching is relatively small, and therefore, negligible. In our research, the exure hinge is modeled as a linear torsional spring rotating about an equivalent rotation 10 axis located at a distance 0:5l h h o from the top end of the rectangular exure element with a constant stiness k h . For the considered wing in passive pitching motion, the deviation rotation is ne- glected. The angular velocity of the wing is a function of apping and pitching angles, ! = _ e x 0 + _ e y = _ cose x + _ e y + _ sine z ; (2.4) where! is expressed in the wing-xed frame. Since the wing is a thin rigid plate, the inertia matrix J w about the xed point o w satises J xz J yz 0 and J z J x +J y . The kinetic energy T and the potential energy V are, respectively, T = 1 2 ! J w ! = 1 2 ! T J w ! = 1 2 ( _ 2 sin 2 + _ 2 )J y + _ _ cosJ xy + 1 2 _ 2 J x ; (2.5) V = 1 2 k h 2 : (2.6) On potential energies, only the hinge elastic energy is considered. The gravitational energy is neglected due to small wing mass and small variation of center of mass position (The gravitational energy is around two orders smaller than the elastic en- ergy). Substituting Eqn. (2.5) into Eqn. (2.1), the equation of motion for the passive pitching of the wing is derived as J y +J xy cos 1 2 J y _ 2 sin 2 = ae @ _ ! @ k h = ae;y k h ; (2.7) where = ae;y k h is the generalized pitching torque about the axisy, including the aerodynamic moment and the hinge restoring torque. Details of the generalized torque calculation can be found in [23]. Employing the same general methodology, the equation of motion for the apping motion is J x +J y sin 2 + _ J y _ sin 2 +J xy cos _ 2 sin = ; (2.8) 11 with = ae;x 0 + drv , where ae;x 0 is the component of ae along axisx 0 of F 0 and drv is the driving torque about the axisx 0 . Since the apping motion is prescribed, the driving torque can be computed from Eqn. (2.8). To determine the attitude of the wing, Eqn. (2.7) is the only equation to be solved. Thus, the pitching motion of the wing is then referred to as the dynamics of a single wing. In Eqn. (2.7), the inertial parameters are estimated from the CAD model, and the hinge is modeled as a torsional spring. To approach the computation of, the only unknown instantaneous aerodynamic torque should be estimated rst. 2.2.2 Integrated dynamics simulation with CFD solver The dynamic equation for the pitching motion of the wing is derived using the La- grangian method in Kane's formulation. To solve Eqn. (2.7), aerodynamic torque must be rst estimated analytically or numerically. During the research initializing procedure, the discretization-based CFD solver was preferred to achieve high-delity simulation. For the low Reynolds ow passing around an insect-sized wing, the gov- erning equations to be discretized are the three dimensional incompressible Navier- Stokes equations [41], r u = 0; @u @t + uru = 1 rp +r 2 u; (2.9) where is the air density, u is the ow velocity, p is the ow pressure, and t is time. Under the weakly coupled FSI assumption, as depicted in Fig. 2.3(a), the CFD solver for Eqn. (2.9) was integrated with the equation of motion (2.7) to predict the pitching angle of the wing. Using the Naiver-Stokes solver, the requirements on the selection of empirical parameters can be minimized compared to the analytical models in Sec. 2.2.3. However, the development of a CFD solver for 3D complex geometry with mov- ing boundaries is not an easy task. To quickly achieve the integrated simulation in 12 Fig. 2.3(a), the nite volume method-based commercial package STAR-CCM+ was employed in the dynamics simulation of a single wing. The algorithm and implemen- tation of our in-house CFD solver for the moving boundary problems are presented in Chapter 4 as the third stage of our research. One challenge in the 3D simulation of apping wings is that the air ow surround- ing the wing is driven by time-varying moving boundaries. To resolve the moving boundaries, the overset grid approach available in STAR-CCM+ is employed. The mesh topology for the single wing simulation is exemplied with the airfoil shown in Fig. 2.3(b). Two mesh blocks and three mesh regions are created. A Cartesian mesh block is generated for the background ow, and a body-tting mesh block is generated (a) (b) (c) (d) (e) Fig. 2.3. Integrated dynamics simulation and overset grids. (a) Data ow for the integrated FSI simulation. (b) Illustration of overset grid assembly with elliptic wing airfoil. (c) Assembly of 3D mesh blocks. (d) Topology for body-tted mesh of the wing with a regular interface. (e) Topology for body-tted mesh of the wing with a non-regular interface. 13 for the wing with high resolution cells near the wall. In general, each mesh block can be assigned either with a prior prescribed or dynamically determined velocity. During the simulation, the two mesh blocks are assembled together dynamically and creating an overlapping region for data exchange. With the overset grid approach, there is no need to re-generate the mesh, and the high quality of the mesh can be maintained in the entire simulation. The accuracy of the overset grid solver can be aected by the interpolation scheme in the overlapping region. The mesh topology for the integrated simulation is depicted in Fig. 2.3(c). The mesh for the background block is not visualized. In the procedure of developing our integrated simulation, three strategies were employed for the generation of the body-tted mesh block: (i) unstructured tetrahedral grid with automated generation algorithm, (ii) structured hexahedral grid with a regular overset interface in Fig. 2.3(d) and (iii) structured hexahedral grid with a non-regular overset interface in Fig. 2.3(e). Strategy (3) provides the highest quality of mesh with the smallest amount of cells. All mesh blocks are generated with ANSYS ICEM CFD. The solver in STAR- CCM+ is set as three-dimensional, laminar, implicit unsteady segregated ow with constant density. Turbulence eects are neglected for the considered low Reynolds number range. The solving process of Eqn. (2.7) is implemented with the transformed set of rst-order equations _ =; _ =F (t;; _ ; ; ae;y ) =F 1 ( ae;y (t)) +F 2 ((t); _ (t); (t)); (2.10) where F 1 = ae;y J 1 y and F 2 =k h J 1 y J xy J 1 y cos + 0:5 _ 2 sin 2. By inspec- tion, it seems that Eqns. (2.10) can be easily solved using standard numerical ODE schemes, such as the 2nd-order modied Euler scheme. However, a major diculty is that the aerodynamic moment ae;y does not have an explicit expression. If a high- order scheme (i.e., higher than the 1st-order scheme) is employed, the expression 14 (a) (b) (c) (d) Fig. 2.4. Quasi-steady components of aerodynamic forces on one wing section for hovering ight. (a) Translational component due to angle of attack and incident velocity. (b) The rst rotational component to satisfy Kutta condition. (c) The second rotational component due to damping. (d) Added mass component due to wing acceleration. of _ cannot be properly evaluated for sub-steps. A compromised approach is the implementation of an incomplete modied Euler scheme, in which the aerodynamic quantity F 1 computed at time t n remains constant from t n to t n+1 . The expression for the intermediate source term is F =F 1 ( ae;y (t n )) +F 2 ( ; _ (t n+1 ); (t n+1 )): (2.11) and F can be thought of as a priori predictions of (t n+1 ) and F (t n+1 ), respec- tively. Using this modication (2.11), the equation set (2.10) is solved in a quasi-2nd order numerical procedure. An extrapolation-based quasi-4th order Runge-Kutta scheme will be discussed in the next chapter. The experimental validation of the integrated simulation is presented in Sec. 2.3.1. 2.2.3 Quasi-steady models of a wing in passive pitching CFD solver can provide us with high-delity estimations for the aerodynamic forces on a wing. However, three drawbacks exist for the integrated simulation with a CFD solver: (1) The pure discrete data from numerical simulation cannot provide us with direct physical insights of geometric and kinematic parameters into the ow phenomena; (2) From a design perspective, a quick estimation in seconds or minutes is expected for iterative design. The integrated simulation for a single wing with a CFD solver typically takes hours or days to compete. As presented in Sec. 3.8, the 15 full dynamics simulation with a CFD solver for a physical ight of 0.2s can take days to nish with 64 CPU cores. (3) We have no understanding of the conditions for the weakly coupled FSI assumption. The integrated simulation with a CFD solver is not computationally ecient for iterative mechanical or controller design. An alternative fast solution must be found. On the modeling of unsteady wing aerodynamic forces, two types can be identied: the quasi-steady contributions as functions of instantaneous kinematics and the highly unsteady contributions that are dependent on the history of kinematics [4]. The quasi-steady contributions dominate the aerodynamic force generation on the wing [6, 7, 20]. Also, there still lack reliable models for highly unsteady contributions. In this study, we primarily employ and adapt existing quasi-steady models for our presented simulations. In Fig. 2.4, four force components are considered to compute the instantaneous aerodynamic force on each wing section in hovering ight [7,20,42{45], df inst = df trans + df rot,1 + df rot,2 + df add ; (2.12) where df trans is the magnitude of the translational force, df rot,1 is the rst subcom- ponent of the rotational force due to circulation, df rot,2 is the second subcomponent of the rotational force due to damping and df add is the magnitude of the added mass force. Here, the viscous shear stress and buoyancy force on the wing surface are neglected [23]. The generation mechanisms for each component are available in [7,20,42,43,45]. The formulas to compute each component are, respectively, df trans = (C L cos +C D sin) 1 2 v 2 n c(r)dr (2.13a) df rot,1 =jv n j rot =C rot jv n j _ c 2 (r)dr; (2.13b) df rot,2 = 1 2 C rd _ _ dr Z te le xjxj dx; (2.13c) df add =C add m ref a a,n = 1 4 C add c 2 (r)a a,n dr; (2.13d) 16 where C L is the lift coecient, C D is the drag coecient, v n is the incident velocity (spanwise ow is neglected), is the local angle of attack (angle between the incident velocity and wing chord),c(r) is the chord length of the wing section at the spanwise location r, rot is the rotational circulation, C rot is the rotational circulation coef- cient, C rd is the rotational damping coecient, C add is the added mass correction coecient, m ref is the reference mass of the section circular disk (Fig. 2.4(d)) and a a,n is the normal acceleration at the midpoint of the wing section. In some cases (e.g., the underwater test shown in Chapter 3), a correction term is needed for f rot,2 to avoid discontinuity of the force curve, the corrected formula for Chapter 3 case is df rot,1 =C rot jv n j _ c 2 (r) sin (2)dr. The correction provides better results than the original form in Eqn. (2.12c). For a 3D wing, the blade element approach is em- ployed to integrate the aerodynamic forces on each wing section along the spanwise direction [28,46]. For a wing moving in the air, the translational force component f trans dominates its aerodynamic force generation. The magnitude of f trans depends on C L , C D , v n and as shown in Fig. 2.4(a). For a single wing in hovering ight, the magnitude of incident velocity and the angle of attack are, respectively, v n =r _ ; (2.14) = 2 sign(v n e x ): (2.15) Here,v n is the characteristic velocity of this wing section and is selected as the velocity at the intersection point between the spanwise axis of the hinge and the wing section. The classic formulas for C L and C D by Dickinson's group are in simple sinusoidal forms [6,7,20,35] as C L =C L 0 sin 2; (2.16a) C D = C Dmax +C D 0 2 C Dmax C D 0 2 cos 2; (2.16b) 17 where C L 0 is the maximum lift coecient, C D 0 is the minimum drag coecient and C Dmax is the maximum drag coecient. These formulas work very well for the cases in which the kinematics of the wing is in a trapezoidal-like shape and the phase dierence between the apping and pitching angles is around =2. The force generation mechanisms for the apping wing at low Reynolds number are complicated. As indicated in [45], the phase dierence between the apping and pitching angles can signicantly alter the leading-edge vortex generation. The variations in the kinematic shape (non-trapezoidal) and dimensionless numbers (e.g., Re) can complicate the situation. Accordingly, the reliability of the lift coecient formula (2.16(a) can decrease. On more sophisticated modeling, [47] incorporated the leading-edge suction analogy into the quasi-steady formulas by adding higher-order sinusoidal terms. [44] employed a parameterized cubic polynomial to represent the asymmetric eects of the lift coecients due to leading vortex variation and other factors. We borrowed the ideas from [44, 47, 48], and added two additional higher- order sinusoidal terms into lift coecient formula of Dickinson's model to enable the capability of the formula to adapt to dierent kinematics. The revised lift coecient formula is C L =C L 1 sin 2 +C L 2 sin 2 cos +C L 3 sin 2 sin (2.17) whereC L 1 corresponds to the primary contribution from the translational velocity and angle of attack,C L 2 andC L 3 re ect unknown unsteady contributions due to dierent kinematic shapes. Based on our simulations results in Sec. 2.3.3, the drag coecient formula is able to provide consistent accuracy for dierent wing kinematics, which is the reason for keeping its original form in our study. On the added mass component, the characteristic acceleration a a,n is the normal acceleration of the disk's center, and the 3D formula proposed in [43] yields a a,n =r cos +x d _ 2 sin cos ; (2.18) 18 wherex d is the distance from the spanwise axis of the hinge to the section's midpoint. Same to [45], C add is introduced as a correction for the cases when the circular disk assumption generates signicant errors. Important aerodynamic moment components include the normal moment inst,n about thex i axis and the pitching moment inst,p about they i axis. Since the viscous shear stress on the wing surface is neglected [23], the moment about the z i axis is assumed to be zero. The moment formulas for the normal and pitching contributions are, respectively, d inst,n = C nc df inst r; (2.19a) d inst;p = d trans,p + d rot,1p + d rot,2p + d add,p = df trans (^ x cp cl h ) + df rot,1p (^ x rot,1 cl h ) 1 2 C rd _ j _ jdr Z te le x 2 jxj dx df add x d ; (2.19b) where C nc is the normal moment correction coecient, d p is the pitching moment corresponding to the section force df , ^ x cp is the normalized location of the center of pressure measured from the leading edge (normalized by local wing chord c), ^ x rot,1 is the acting location of the rst rotational force, and x d is the distance from the hinge axis to the section midpoint. Calculations of d inst,n , d trans,p and d rot,2p are straightforward with variables and parameters from Eqn. (2.13). d trans,p and d rot,1p require models for the dimensionless parameters ^ x cp and ^ x rot,1 . For ^ x cp , as discussed in [7,28,32,47,49], broadly accepted models include the quarter-chord model ^ x cp = 0:25 and the empirical model ^ x cp = 0:82 jj + 0:05: (2.20) The above empirical model is referred to as the \classic model" later. The quarter- chord model comes from the thin airfoil theory (small angle of attack and high 19 Reynolds number). The classic model considers the variation of angle of attack and can provide more information about the movement of the center of pressure. However, the experiments for the empirical model were conducted statically and the angle of attack was unchanged for the region to characterize the measured aerodynamic force. As shown in [47] (see its Fig. 10), the static empirical model deviates signicantly from the results of dynamic experiments in which the pitching angular velocity is a nonzero constant. From Fig. 10 in [47], we can also see that the trend line of ^ x cp for angles between =6 and =2 passes through the point (0, 0:25c), which is consistent with the quarter-chord model. Without introducing new mechanism for the rota- tional force [7], we employ a model consistent with the dynamic experiments in [47] by adjusting the slope and intercept of the classic formula ^ x cp = 0:4 jj + 0:25: (2.21) For = =2, ^ x cp is 0:45 not 0.5 because the wing planform is not rectangular. The coecient 0:4= can be adjusted according to dierent wing planforms. Notice that the center of pressure in [47] accounts for contributions from all force components, but ^ x cp in Eqn. (2.21) only considers the contribution from the translational force. According to [47], the midpoint model of ^ x rot,1 is used, i.e. ^ x rot,1 = 0:5. To summarize, the coecients C L 1 , C L 2 , C L 3 , C Dmax , C D 0 , C rot , C rd , C add and C nc form the dimensionless parameter set to calibrate for the accuracy adjustment of the aerodynamic force calculation. These dimensionless parameters are calibrated using experimental or CFD data for a 3D wing, and the initial values can be chosen from the classical results in [7,20,35,47]. 20 d l 0 1 (a) (b) (c) β θ θ limit θ limit Fig. 2.5. Geometric limit by hinge shoulders and correction methods. (a) Geometric limit applied by hinge shoulders. (b) Geometric limiter method. (c) Enlarged stiness and damping method. 2.2.4 Pitching angle limit by hinge shoulders and wing de- formation In the practical design of a exure hinge, there is only a limited number of thickness values available for the Kapton lm. For dierent wing shapes and structures, the enter of mass location and moment of inertia distribution are all dierent. Accord- ingly, to constrain the maximum pitching angle of the wing, it is not always feasible to avoid a geometric collision illustrated in Fig. 2.5(a). If the collision happens, there are two physical phenomena: (1) The carbon ber presents small deformation to permit a continuing increase of the pitching angle. (2) The deformation of the Kapton lm changes into tension mode from bending mode. Both cases provide a signicantly enlarged stiness. Since the pitching angular velocity is not zero when the collision happens, the angular velocity is damped by the collision quickly. Since the experimental investigation of the physical procedure for the collision is not easy, we only provide some numerical strategies to mimic the behavior of the collision pro- cess. Two methods were tested in our numerical simulations: the geometric limiter method and the enlarged stiness-damping method. In the geometric limiter method, the numerically updated pitching angular veloc- 21 ity with Eqn. (2.7) is corrected by a limiter coecient for each time step size, _ = _ comp ; (2.22) where _ comp is the computed value with Eqn. (2.7), and the limiter is described in Fig. 2.5(b). The limit value by the shoulder geometry is limit limit and limit re ects the rigidity increase by the collision. The results in Sec. 2.3.5 will show the disadvantage of this sharp correction prole and smoother prole is expected in practical use of this method. In the enlarged stiness-damping method as in Fig. 2.5(c), extra torque term re ecting the angle limit is added into Eqn. (2.7), extra =k c ( limit ) +c c _ ; > limit > 0; (2.23) where k c and c c are the correction stiness and damping. Eqn. (2.23) is > 0, a negative sign is applied to the right-handed side of Eqn. (2.23) for < 0. Since the correction stiness and damping cannot be experimentally obtained in a simple method, the values in our simulations are estimated based on the observations from the maximum pitching angle in the tests. In [22, 24, 50{53], the wing deformation is also modeled as the passive rotation with a stiness. In some cases, such as underwater tests in Sec. 3.8.4, even the pitching angle of the wing is limited by the geometric dimensions, the wing can still pitch further with deformation. Therefore, Eqn. (2.23) is also the simple model for wing deformation. Dierent from the collision case, the correction stiness for wing deformation can be a negative value. The overall estimated values of k c and c c for collision and wing deformation are used in Sec. 3.8.4. 22 2.2.5 Amended model for the dynamics of a wing in water In the next chapter, open-loop experimental tests in the air were run for the valida- tion of our full dynamics model, but ended in unsatisfactory results from quantitative purpose. We turned into the underwater tests described in Sec. 3.6. However, if the process in Sec. 2.1 2.3 is directly applied to the underwater tests, no converged simulations were achieved for either CFD or QS models. Multiple corrections, in- cluding stier hinge, larger wing moment of inertia, and dierent ODE schemes, were tried but failed without little help. A quantitative analysis was made using the QS models to identify the reason. The problem arises from the added mass term in the lumped aerodynamic torque in Eqn. (2.7). We decompose the lumped aerodynamic torque into two parts: ae;y = ae;y=add J add ; (2.24) whereJ add is the equivalent pitching moment of inertia from added mass force. In our experimental tests, the apping frequency for the case in the air is around 100 Hz, and the frequency for the case in the water is around 13 Hz. The pitching moment of inertia of the wing J y is around 3 mg mm 2 , the value of J add in the air at 100 Hz is around 0:8 mg mm 2 , and the value ofJ add in water at 13 Hz is around 640 mg mm 2 . From the perspective of numerical scheme, the correction value of 0:8 cannot change the convergence property of Eqn. (2.7), but the correction value of 640 is able to scale-up any numerical errors very quickly. The pitching motion equation (2.7) with the lumped aerodynamic torque must be modied to guarantee converged numerical results. The amended model with J add is (J y +J add ) +J xy cos 1 2 J y _ 2 sin 2 = ae;y=add k h ; (2.25) where J y +J add is the augmented moment of inertia for the wing. Since J add is 23 (a) Measured angle trajectories (b) (c) Time(ms) 0 5 10 15 20 0 50 -50 R = 15 mm Fig. 2.6. Experimental setup and measured kinematics from [7] for validation of the single wing dynamics. (a) Illustration of the experimental setup. (b) Planar shape of the wing in experiment and simulations. (c) Experimentally measured apping angle , pitching angle and deviation angle in [7]. signicantly larger than J y , it is reasonable that the adjustment of J y does not help for the convergence. The converged results in Sec. 2.3.6 demonstrated our analysis and eectiveness of our amended. In the integrated simulation with CFD solver, the added mass is not able to be directly separated from the simulation. To enable a converged CFD simulation, an estimated value from the QS model can be employed because the exact value of J add is not important. The amended form presented in (2.25) was rst derived in [54] to calculate the resonant frequency of the pitching motion. But, the importance of the relative mag- nitude between J y and J add was not pointed out for simulation in a denser media. The amended model extends the adaptability of our integrated simulation. More in- terestingly, this provides us with an example for the validity of the weakly coupled FSI assumption. 2.3 Results 2.3.1 Experimental validation The integrated dynamics simulation with a CFD solver described in Fig. 2.3(a) was validated by comparing the computed pitching angles and instantaneous vertical force with the experimental data available in [7]. The benchmark results from [7] 24 -2 0 2 0 20 15 10 5 (b) Inertial force components in vertical direction Time (ms) Force (mN) -2 0 2 (c) Total vertical force with Fourier series fitting Time (ms) Force (mN) 0 20 15 10 5 (d) -2 0 2 Total vertical force with zero-phase filter Time (ms) 0 20 15 10 5 Force (mN) Angle (o) (a) 0 15 10 5 30 25 -50 50 Simulated and experimental pitching angles Time (ms) 0 20 Fig. 2.7. Model validation for the dynamics of passive wing pitching at 100 Hz. (a) Instantaneous pitching angles comparison. (b) Vertical inertial force computed from simulated and measured . (c)-(d) Total vertical force with dierent ltering methods for measured . FS represents Fourier series tting and ZP denotes zero- phase ltering. were obtained through the experimental setup graphically described in Fig. 2.6(a), which consists of the object wing, custom-made force sensors, high-speed cameras, and optical components. The articial wing is biologically inspired by the species Eristalis tenax (drone y). The planar shape of the wing is formed by the polyimide lm, and the wing is structurally reinforced with carbon ber as veins to make it as rigid as possible. From the green line in Fig. 2.6(c), the deviation angle is at least one order smaller than the apping (blue) and pitching (red) angles, which is consistent with our assumption that _ 0. However, as demonstrated in [23], the variation of in high frequency can signicantly aect the estimation of vertical inertial force. The experimentally measured deviation angles must be included when computing the vertical inertial force. In the benchmark experiment and simulations, the wing aps in a sinusoidal-like prole, which is not perfectly symmetric, as shown in Fig. 2.6(c). The initial con- ditions and inertial properties of the wing are listed in Table 2.1. The simulated 25 pitching angles present an excellent agreement with the experimental ones (in terms of magnitudes and phases) in most regions except for the peaks marked with A. The angle matching directly demonstrates the eectiveness of the proposed integrated sim- ulation approach. In Fig. 2.7(b), the vertical inertial force components are calculated from wing kinematics using formulas presented in [23]. The pitching-dominated iner- tial force ma is computed with the simulated pitching angles (blue curve in 2.7(a)) while the deviation dominated inertial force ma is estimated with the experimental deviation angles (green curve in 2.6c)). Since the measured deviation angles are in small magnitudes with high-frequency oscillations, the experimental data is processed with two approaches: Fourier series tting and zero-phase ltering. Due to the high- frequency oscillations, the inertial force from the deviation angle is in the same order to that from the pitching angle. The deviation angle cannot be neglected in total force calculation. The resulted total vertical forces are presented in Figs. 2.7(c)-(d). Both processing strategies are able to provide good agreement with the instantaneous experimental force. However, in regions C and D, the two strategies provide dierent levels of agreement. Since the deviation angle is reconstructed from images instead of being directly measured, the discrepancies on the instantaneous matching (e.g., region E) could be attributed to the measurement errors of the deviation angle. From the mean value perspective, the measured cycle-averaged vertical force f L in [7] is 71.6 mg. The simulated cycle-averaged aerodynamic force f L;ae is 69.1 mg. Theoretically, the average inertial force in one cycle should be zero for a periodic motion. The simulated cycle-averaged vertical force f L is 69.1 mg. The relative error between the simulated and experimental results is 3:5%. Considering that the custom-made force Table 2.1: Parameters of the wing for the benchmark passive simulation. (0) (0) k h t J x J y J z J xy -16.17 54.85 2.352 0.1 43.0 1.7 41.3 -3.5 Note: The units for angle, stiness, time and moment of inertia are , Nm, ms and mg mm 2 , respectively. 26 sensor is very sensitive to the operating conditions [7], this level of accuracy is a very satisfactory agreement. In the validation procedure, we did not explicitly tune any parameters. It indicates the dynamics model can truly re ect the physics of the wing pitching. In Fig. 2.7(a), a noticeable discrepancy between the simulated and experimen- tal pitching angles is observed in region A, where the simulated angle has a larger magnitude and also exhibits a slight phase lag. Possible reasons for the discrepancy include the computing error of the aerodynamic force, the nonlinear behaviors of the exure hinge, the measurement error of the pitching angle, etc. On the force errors, a direct measurement of the instantaneous aerodynamic force instead of the total vertical force is required. On the hinge modeling, as discussed in Sec. 2.2.1, the linear model for the hinge stiness may not work for a large pitching angle because of geometric nonlinearity. Also, there possibly exist geometric collisions between the shoulders of the hinge to limit the pitching angle as discussions in Sec. 2.2.4. About the measurement errors, since the camera is put at one side of the wing instead of from the top, the near and far positions can generate dierent levels of errors for the pitching angle estimation. More precise experiments for wing kinematics measure- ments are required to validate our guesses further. In summary, despite the minor short-duration discrepancies, the integrated simulation can provide us with accurate enough results for the design purpose. 2.3.2 Applications of the integrated simulation with CFD With the proposed integrated dynamics simulation, we can investigate the aerody- namic force generation on the wing with dierent wing kinematics, geometry and inertia. As an example, we simulated the same wing with various values of hinge stiness. The calculated pitching angles and vertical aerodynamic forces are pre- sented in Figs. 2.8(a)-(b). When the hinge stiness decreases from 4 to 1 Nm, the amplitude of the pitching angle increases and the peak shifts right. Accordingly, the 27 (a) 0 -50 50 θ for various hinge stiffnes 0 5 20 15 10 25 Times (ms) θ (o) Times (ms) 0 2 f L,ae for various hinge stiffness 0 5 20 15 10 25 f L (mN) (b) Honeybee Eristalis tenax θ for different wing shapes 0 5 20 15 10 25 Times (ms) θ (o) 0 50 -50 100 (c) Honeybee Eristalis tenax f L,ae for different wing shapes 0 5 20 15 10 25 Times (ms) 0 2 f L (mN) (d) Fig. 2.8. Parametric studies of hinge stiness and wing planar shape. (a)-(b) Pitch- ing angles and vertical aerodynamic forces for various hinge stiness. (c)-(d) Pitching angles and vertical aerodynamic forces for a Eristalis-tenax-like wing and a honeybee- like wing for a hinge stiness k h = 3. The unit for hinge stiness is Nm. Table 2.2: Cycle-averaged lift for various hinge stiness. k h (Nm) F L,ae (N=mg), 100Hz F L,ae (N/mg), 120Hz 1 363.6/37.1 216.6/22.1 2 669.3/68.3 715.4/73.0 3 805.6/82.2 968.2/98.8 4 837.9/85.5 1005.6/102.6 maximum value and amplitude of the vertical aerodynamic force become smaller as k h decreases. For current wing shape and apping prole, the hinge stiness 4 Nm is the best choice in terms of lift force generation. The cycle-averaged lift for two frequencies are listed in Table 2.2. The results for the cycle-averaged lift are consis- tent with the instantaneous curves. After calculating the growth rate of f L,ae for each frequency, the case at 120 Hz is more sensitive to the value of hinge stiness. The simulation results for two dierent wing planar shapes are presented in Fig. 2.8(c)-(d), one is a Eristalis-tenax-like wing and the other one is a honeybee-like wing. 28 The shapes of the two wings are adjusted so that the two wings have very similar spanwise lengths and surface areas. The moment of inertia is also assumed the same for both wings, which can be achieved by adding lumped masses onto the wings. The main dierence is that the honeybee-like wing has a sharper trailing edge than the Eristalis-tenax-like wing. From Fig. 2.8(c)-(d), compared to the other wing, the honeybee-like wing has a signicant phase shift in the pitching angle prole, and larger instantaneous vertical forces, which proves that the planar shape can signicantly aect the aerodynamic eciency of the passively pitching wing. Also, the moment of inertia of the wing can be further adjusted to achieve better aerodynamic performance. 2.3.3 Adaptability of the revised QS models The adaptability of our revised quasi-steady models for a single wing was numerically demonstrated with CFD simulations. To precisely describe the kinematics of the wing, the active simulations with fully prescribed wing kinematics are employed as = 0 sin(2t) and = 0 sin(2t + t ). For the results in Fig. 2.9, the simulations were implemented at two t values (0.5 and ) with = 100 Hz, 0 = =3 and 0 = 0:995. Three sets of parameters (i.e., P 1 , P 2 and P 3 ) are tested to illustrate the reliability of the adapted models. The values of these parameters are listed in Table 2.3. These parameter sets are calibrated based on the data agreements for the lift and drag curves rather than the moment curve. P 1 strictly satises Eqn. (2.16)(a) while P 2 follows Eqn. (2.17) with a non-zero C L 2 . In Fig. 2.9(a), the dierences between P 1 and P 2 are not signicant for the Table 2.3: Parameter sets for the quasi-steady models of a single wing. Name C L 1 C L 2 C L 3 C Dmax C D 0 C rot C rd C add C nc ^ x cp P 1 1.35 0.00 0.0 2.38 0.29 2.5 2.0 1.75 0.9 (2.21) P 2 1.17 0.25 0.0 2.38 0.29 2.5 2.0 1.75 0.9 (2.21) P 3 1.17 0.25 0.0 2.38 0.29 2.5 2.0 1.75 0.9 (2.20) 29 QS-P 3 CFD QS-P 1 QS-P 2 (a) -4 0 4 Pitching moment for Δθ t = π -1 0 1 Pitching moment for Δθ t = 0.5π Time (ms) 0 10 20 30 40 Time (ms) 0 10 20 30 40 Time (ms) 0 10 20 30 40 0 1 2 Lift for Δθ t = 0.5π f L (mN) (b) -2 0 2 Drag for Δθ t = 0.5π Time (ms) 0 10 20 30 40 f D (mN) τ p (μNm) τ p (μNm) (c) -1 0 1 2 Lift for Δθ t = π Time (ms) 0 10 20 30 40 f L (mN) (d) -4 0 4 Drag for Δθ t = π Time (ms) 0 10 20 30 40 f D (mN) (e) (f) Fig. 2.9. CFD validation of quasi-steady models for a single wing in hovering. (a)- (c) Instantaneous lift, drag and pitching torque for t = 0:5. (d)-(f) Instantaneous lift, drag and pitching torque for t =. t = 0:5 case. However, in Fig. 2.9(e), P 1 has larger discrepancy than P 2 at the low lift regions for the t = case. The amplitude error is around 20%. More importantly, the larger discrepancy with P 1 cannot be removed by only tuning C L 1 because the discrepancy in other regions can then increase. The drag coecient formula (2.16)(b) provides consistent accuracy for both t cases in Fig. 2.9(b) and (d), which demonstrates the adaptability of the classic drag formula. The results for parameter sets P 2 and P 3 are shown in Fig. 2.9(c) and (f). Notice that the parameters of P 2 calibrated with the lift and drag curves automatically satisfy 30 (a) Pitching angle Lift Pitching moment Drag 0 10 20 30 40 50 60 Time (ms) 0 10 20 30 40 50 60 Time (ms) 0 10 20 30 40 50 60 Time (ms) 0 10 20 30 40 50 60 Time (ms) θ (rad) 2 -2 0 τ p (μNm) 2 -2 0 f L (mN) 0 2 f D (mN) 4 0 -4 (b) (c) (d) Fig. 2.10. Quasi-steady model calibration for a single wing in passive pitching. (a)-(d) are instantaneous pitching angle, pitching moment, lift and drag, respectively. the accuracy of the pitching moment estimation for both t cases. This is a strong evidence for the adaptability of the revised models (2.17) and (2.21). In contrast, P 3 , which follows Eqn. (2.20), does not correctly estimate the pitching moment for t = 0:5 in Fig. 2.9(c). P 3 results in smaller pitching moment than the CFD data in most regions. The above dierence between P 2 and P 3 is consistent with the under-damped pitching behavior observed in [7], in which the model P 3 is applied. The above active simulations demonstrate our revised models for the translational force component present improved adaptability for dierent wing kinematics, which is crucial for the full dynamics investigation of the entire MAV using QS models in the next chapter. 2.3.4 Integrated dynamics simulation with QS models For a wing in passive pitching, the integrated dynamics simulation was carried out for the calibration of the quasi-steady models in Sec. 2.2.3 with the results in Fig. 31 Table 2.4: Parameter sets for the quasi-steady model of a wing in passive pitching. Name C L 1 C L 2 C L 3 C Dmax C D 0 C rot C rd C add C nc ^ x cp k h l h P 4 1.00 0.85 0.00 3.40 0.40 1.0 2.0 1.0 0.9 (2.21) 1.5 0.8 P 5 0.85 0.72 0.00 2.72 0.32 1.6 2.0 1.0 0.9 (2.21) 1.5 0.8 Note: The unit for k h is Nm and the unit for l h is mm. 2.10. The parameter set for the quasi-steady model is P 4 in Table 2.4. The apping prole of the wing is(t) = 0 sin(2t) with 0 ==3. As displayed in Fig. 2.10(a), after the ow reaches a periodic steady state, the pitching angles predicted by the quasi-steady model are almost identical to those of the CFD simulations. According to Eqn. (2.7), when the inertial properties, hinge stiness and apping prole of the wing are given, the pitching angle is determined only by the aerodynamic moment. Thus the agreement of the pitching moment is consistent with the results of the pitching angle. Similar to the pitching moment curves, the drag curves also have exceptionally satisfactory matching to each other. The lift curves show good but imperfect instantaneous matching. The largest discrepancies exist near the peak and valley regions. However, if the data are averaged from the 3rd to 10th period, the average lift of the CFD simulation is 63.5 mg, and the average lift of the quasi-steady analysis is 64.1 mg. The relative error is less than 1%. The cycle-averaged lift is in high accuracy. Since the instantaneous lift estimation is good in most regions and the matching for instantaneous moment and drag are satisfactory, the errors of the estimation of the forces shown in Fig. 2.10 are assumed acceptable for our study in this session. In Sec. 2.2.3, there are a set of parameters that can be tuned for the calibration of the QS models of a wing. Since there is not a unique standard to select those parameters, it is highly possible to nd more than one set of parameters, e.g., P5 in Table 2.4. For the lift curves, the matching in the extreme regions can be improved with P5 in which the coecient C rot is modied to a higher value. The parameter set P5 is able to provide a similar order of accuracy for the given wing and apping 32 prole. However, we cannot directly tell which set of parameters provides a better representation of the physics for force generation on the wing, or which set of param- eters is only for purely data tting. A more robust calibration method is discussed in the full dynamics simulation in Sec. 3.8.4 2.3.5 Pitching angle limit by hinge shoulders Pitching angle comparison Time (s) 0 0.1 0.2 0.3 0.4 -1 0 1 θ (rad) -2 0 2 (rad/s) 10 5 Pitching angluar acceleration Time (s) 0 0.1 0.2 0.3 0.4 (a) (b) limiter free Fig. 2.11. Pitching angle limit with the geometric limiter method. The results with the quasi-steady models and geometric limiter method are plotted in Fig. 2.11. The parameters for the permitted pitching angle are limit = =6 and = 0:01. From Fig. 2.11(a), the limiter method does eectively constrain the maximum pitching angle of the wing. However, as seen in Fig. 2.11(b), when the collision happens, there are spikes in the angular acceleration curves. As these spikes can make the added mass force super larger, non-physical behavior can be observed for the full dynamic simulation in the next chapter. Thus, to guarantee the parameters physically reasonable, limit and limit must be calibrated with some experimental tests. The results of the enlarged stiness-damping method resemble the curves in Fig. 2.11. Since the geometric limiter method enforces a hard limit on the pitching angle, it cannot naturally include the wing deformation into consideration. The enlarged stiness-damping method was preferred in our underwater tests in Sec. 3.6. 33 2.3.6 Integrated simulation for a wing immersed under water (a) Pitching angle comparison Time (s) θ (rad) 0 0.1 0.2 0.3 0.4 -1 0 1 θ quasi θ cfd (b) Fig. 2.12. Converged integrated simulations for underwater tests at 13 Hz. (a) Wing planar shape. (b) Converged pitching angle proles. The integrated simulations with both CFD solver and QS models were imple- mented with the amended dynamic model (2.25) and the planar wing shape in Fig. 2.12(a). The apping prole is = 0 sin(2t) with 0 = 0:4436 and = 13 Hz. With the correction moment of inertia J add estimated from the quasi-steady model, the pitching motion solved with CFD solver is able to converge into periodic oscil- lations. The CFD and QS models provide satisfactory matching with each other on the prediction of the pitching angle, which is consistent with the results in Sec. 2.3.4 for the case in the air. 34 Chapter 3 Full Dynamics of the Entire Flapping-Wing MAV 3.1 Problem description The dynamics model of a single wing provides us with an eective tool to analyze the aerodynamic eciency and select geometric or inertial parameters for the wing. However, to answer questions such as guidelines for the installation location of the wing onto the body and the eects of the wing inertia on the initial response of the MAV, dynamics model for the entire MAV must be formulated. For the entire MAV, its free ight in space is also a FSI process as shown in Fig. 2.1(c). The air ow generates aerodynamic forces on the body and wings to change the motion of the MAV while the movement of the body and wings act as moving boundaries to drive the ow eld to evolve. The dynamic behavior of the MAV is the result of the interactions between the inertial, elastic and aerodynamic forces. To avoid numerical diculties and make the integrated analysis easier to implement, the weakly decoupled FSI assumption is also applied to the entire MAV analysis. In this chapter, we are to derive the equations of motion for the entire MAV to achieve the capability of the simulating the free ight of the robotic MAV as depicted in Fig. 3.1(a). To provide kinematic descriptions for the robotic MAV, six right-handed frames 35 Trajectory MA V (a) Left wing Right wing Top view Left side view (c) (b) Fig. 3.1. Free ight and kinematic description of the MAV. (a) An illustration of the free ight. (b) Frames of reference. (c) Flapping and pitching angles with respect to the body. of reference are employed as shown in Fig. 3.1(b) (Two frames of reference for the other wing are not plotted for gure simplicity). The inertial frameF 0 fo 0 ;x 0 ;y 0 ;z 0 g has its origino 0 specied as the MAV's initial position, itsx 0 axis pointing vertically upward and itsy 0 axis along a specied horizontal direction. The body-xed frame F b fo b ;x b ;y b ;z b g denes its origino b as the center of mass of the MAV body, itsx b axis pointing along its body axis toward the MAV's head and itsy b axis normal to the body axis and pointing rightward. The frame F r;i fo r;i ;x r;i ;y r;i ;z r;i g is located at the root of the wing i (i = 1; 2). Thex r;i axis is normal to the stroke plane of the wing andy r;i is normal to the symmetric plane of the body. In the presented design, the wing stroke plane is always normal to the body axis, so the three axes of F r;i are parallel to the three axes of F b . Thus, F r;i is also a body-xed frame. The wing-xed frame F i fo i ;x i ;y i ;z i g is located at the center of mass of the wing i with itsy i axis parallel to the spanwise axis of the hinge and pointing toward rightward, itsz i axis normal to the wing plane and itsx i pointing to the leading edge of the wing i. The initial orientations for all the frames of reference are identical. Translational motion of the MAV is dened by the relative position of o b with respect to o 0 . The attitude of the MAV is depicted using three consecutive Euler angles between F b and F 0 , that is, roll angle b aboutz 0 0 (parallel toz 0 and located at the center of mass of the body) axis, pitch angle b about the carried axisy 0 0 (y 0 36 after rotating about z 0 0 ) and yaw angle b about the x b axis. The attitude of the wing with respect to the body is depicted using the relative Euler angles between F i and F r;i , that is, the apping angle i about thex r;i axis, the deviation angle i about carried axisz 0 r;i (z r;i after the rotating about x r;i ) and the pitching angle i about the y i axis. Accordingly, the MAV has three translational DOFs and seven rotational DOFs. Two sets of rotations are dened: roll, pitch, yaw for the body and apping, deviation, and pitching for the wing. For the design in Fig. 3.1(a), there is no specic mechanism for the deviation motion except small elastic eects, thus the deviation motion is assumed to be zero for simplicity [23]. Denitions of the maximum, minimum, and zero values for the apping and pitching angles with respect to the body are depicted in Fig. 3.1(c). The concept of \symmetric apping prole" for a sinusoidal signal is dened as the mean position lies inside thex r;i y r;i plane without bias. From the perspective of dynamics simulation, under some conditions, manually deriving equations of motion become not necessary nowadays because the automated symbolic derivation tools (e.g. Maple) or general-purpose solvers (e.g. MSC.Adams) are possibly available for desired studies. However, three reasons motivated us to manually formulate the equations of motion for the entire MAV. First, if the apping motion of the wing is at high frequency, the studied MAV would be a multiple-time- scale system. Small time step size is required for a reliable computation of the instan- taneous wing kinematics. The low computing eciency of symbolic operations can prevent its applicability to fast iterative design or long-endurance ight simulations. Second, the objectives of our dynamics study are not limited to the steady-state and transient response analyses. Physical insights and simplied models are needed for theoretical analysis and parametric design as demonstrated in [27]. The automated derivation or a general-purpose solver cannot provide necessary details on the struc- ture of the dynamic equations. Third, the dynamics simulation for MAVs requires calculation of aerodynamic forces, but it is dicult to strongly couple proprietary 37 dynamics software with aerodynamic solvers, especially CFD solvers. The manually derived equations of motion outperform the automated tools in above three aspects. However, the manual derivation of the equations of motion for the entire MAV is not a straightforward process like that of a single wing. The cumbersome algebraic work can prevent us from obtaining a concise set of equations to conduct theoretical analysis and implement fast numerical computation. Therefore, in this chapter, we rst describe our proposed alternative form of the Lagrange's equations so that a computationally ecient form of the equations of motion can be manually derived for the entire MAV. 3.2 Lagrange's equations alternative The amount of cumbersome algebraic operations prevents Lagrange's equations from being used as an ecient approach to manually derive equations of motion for a complex system [29, 55]. This fact is partially due to the repeated terms existing in the left-hand side of Eqn. (2.1) [39]. During the derivation of Lagrange's equations, these repeated terms are introduced to present a more compact and generalized form of the dynamic equations. If the generalized coordinates are properly selected, these repeated terms can be easily removed. To describe spatial vectors in various reference frames, we usee to represent the unit vector of an axis and e to represent a column set (vector) of the base vectors of a frame F . For example, 0 e = [e x 0 ;e y 0 ;e z 0 ] T is a basis set for F 0 . We can express the position vectorr i = 0 e T 0 r i where 0 r i is the coordinate vector ofr i in frame F 0 . Ifv i and! i are expressed in F 0 and F i respectively,v i = 0 e T 0 v i and! i = i e T i ! i . If J i is the moment of inertia tensor of the ith component, the coordinate matrix of J i inF i is i J i . For simplicity,r i ,v i ,! i and J i are the shorthands for 0 r i , 0 v i , b ! i and i J i respectively. For the discrete system in consideration, the kinetic energy is T = Nc X i=1 1 2 m i v i v i + 1 2 ! i J i ! i = Nc X i=1 1 2 v T i m i v i + 1 2 ! T i J i ! i : (3.1) 38 Here T is a scalar function which depends on the generalized coordinate vector q = [q 1 ;q 2 ;:::;q Nq ] T , the generalized velocity vector _ q and timet, that is,T =T (q; _ q;t). Using the matrix dierentiation conventions dened in [56] and the symmetric prop- erty of J i , we write @T @q = Nc X i=1 v T i m i @v i @q +! T i J i @! i @q (3.2) and d dt @T @ _ q = Nc X i=1 dv i dt T m i @v i @ _ q +v T i m i d dt @v i @ _ q + d! i dt T J i @! i @ _ q +! T i J i d dt @! i @ _ q : (3.3) Since! i might be contained in the expression of velocityv i , the termv T i m i @v i @q in Eqn. (3.2) and the term v T i m i d dt @v i @ _ q in Eqn. (3.3) are highly complicated compared to the other terms. If these two terms are removed, the required algebraic labor can be signicantly reduced in the formulation of the generalized inertial force d dt @T @ _ q @T @q . Fortunately, sincev i is the coordinate ofv i in the inertial frameF 0 , we can have that @v i @q and d dt @v i @ _ q are equal to each other. For the system considered, we assume r i ,v i andq satisfy required continuity and smoothness conditions. Elements ofq are independent of each other as dened in Eqn. (2.1). As r i is the position vector, r i is the function of the generalized coordinates and time, not dependent on the generalized velocities, that is, r i =r i (q;t). Since r i and v i are coordinate vectors in the inertial frame, then v i = dr i dt = @r i @t + s X j=1 @r i @q j _ q j : (3.4) Following this, we nd that @v i @q j = @ @q j @r i @t + s X k=1 @ @q j @r i @q k _ q k = @ @t @r i @q j + s X k=1 @ @q k @r i @q j _ q k = d dt @r i @q j = d dt @v i @ _ q j : (3.5) 39 Note that the equality @r i @q j = @v i @ _ q j is obtained by taking partial derivatives on both sides of Eqn. (3.4) with respect to _ q j [57]. Accordingly, v T i m i @v i @q =v T i m i d dt @v i @ _ q : (3.6) This equality does not necessarily hold for i v i because similar condition presented in Eqn. (3.4) may not exist. Therefore, at least for v i and ! i , the generalized inertial force can be simplied as d dt @T @ _ q @T @q T = Nc X i=1 @v i @ _ q T m i dv i dt + @! i @ _ q T J i d! i dt + d dt @! i @ _ q @! i @q T J i ! i : (3.7) The transpose operation is applied in Eqn. (3.7) because it is easier to obtain the state space representation. The canceling out process is a reversal procedure of deriving Lagrange's equations in [58]. If the reference point for kinetic energy T is not at the center of mass of each component [39], the cancellation property still exists. But there exists coupling terms between the translational velocity and the angular velocity in Eqn. (3.1), which requires additional processing for derivatives. Here, we only illustrate the formulation of dynamic equations with reference points selected at center of mass. Typically, potential energy V is a function of q and t. Because the potential energy only exists in the second term of Eqn. (2.1), there is no similar cancellation available for the potential force. The original form is retained as @V @q T . For the non- conservative forces, a more general formula than the one in Eqn. (2.1) can be obtained because the reference point for the resultant force and moment is not necessary at the center of mass of each component. If P i is a xed point on the ith component, the non-conservative force vector is [23,58] p q = Nc X i=1 @v P i @ _ q T f i + @! i @ _ q T P i : (3.8) 40 Here, f i is the coordinate of f i in F 0 , and P i is the coordinate of P i in F i . If we integrate Eqn. (3.8) with (3.7) and @V @q T , the Lagrange's equations alternative is obtained. This alternative can be treated as an intermediate result when deriving Lagrange's equations. 3.3 Equations of motion The alternative form of the Lagrange's equations is derived in the last subsection, the basic task here is to nd kinematic expressions and corresponding derivatives for each component of the MAV system. The studied MAV has a total of ten DOFs for the body and wings. The position and translational velocity of the center of mass of the body are r b = 0 e T r b ;v b = 0 e T v b = 0 e T _ r b (3.9) where r b = [x b ;y b ;z b ] T is the coordinate vector ofr i in F 0 and v b is the coordinate of v b in F 0 . x b , y b , z b are the three generalized coordinates for the translational motion of the body. Using the Euler angles dened for the body rotation, the angular velocity of the body is ! b = _ b e z 0 0 + _ b e y 0 0 + _ b e x b = b e T 2 6 6 6 4 _ b sin b _ b cos b sin b _ b + cos b _ b cos b cos b _ b sin b _ b 3 7 7 7 5 = b e T T b _ b = b e T ! b ; (3.10) where b e is the basis set ofF b and! b is the coordinate vector of! b inF b . We can see that! b = T b _ b . The Euler angles b , b and b are the three generalized coordinates for the rotation of the body. The transformation matrix T b and the coordinate vector 41 are introduced as auxiliary matrices T b = 2 6 6 6 4 1 0 sin b 0 cos b cos b sin b 0 sin b cos b cos b 3 7 7 7 5 ; b = 2 6 6 6 4 b b b 3 7 7 7 5 : (3.11) Before nding expressions for the wing kinematics, rotation matrices must be dened to convert coordinates between various reference frames. Rotation of the body with respect to the inertial frame can be represented by the relation between b e and 0 e as b e = R b ( b ; b ; b ) 0 e = R b 0 e; (3.12) where R b is the rotation matrix that transforms coordinates between F 0 and F i . Similarly, we can dene the rotation matrix R i as i e = R i ( i ; i ; i ) b e = R i b e; (3.13) where i e is the basis set of F i . For the current case, i is always zero, so i and i are chosen to be the generalized coordinates that describe the relative rotation of the ith wing with respect to the body. Entries of R b and R i are listed in Appendix A. The rotation of the wing with respect to the inertial frame can be decomposed into two sequential rotations, that is, the rotation along with the body and the relative rotation with respect to the body. Using the angular velocity addition theorem [59], the angular velocity of the ith wing is ! i =! b +! i;b =! b + _ i e x r;i + _ i e y i = b e T T b _ b + i e T T i _ i = i e T R i T b _ b + T i _ i = i e T ! i (3.14) where ! i is the coordinate vector of ! i in F i . The transformation matrix T i and 42 coordinate vector i are dened as T i = 2 6 6 6 4 cos i 0 0 1 sin i 0 3 7 7 7 5 ; i = 2 4 i i 3 5 : (3.15) The above matrices are for the 2DOF wing rotation case and they can be easily extended into the 3DOF case when non-zero i is considered. As shown in Fig. 2.1, the position vector of the center of mass of the ith wing is r i =r b +r r;i +r c;i = 0 e T r b + b e T r r;i + i e T r c;i = 0 e T r b + R T b r r;i + R T b R T i r c;i = 0 e T r i : (3.16) wherer b is the coordinate ofr b inF 0 ,r r;i is the coordinate ofr r;i inF b andr c;i is the coordinate ofr c;i in F i . Since the body and wings are rigid, r r;i and r c;i are constant coordinate vectors. The velocity vector of the center of mass of the ith component is v i = _ r i = 0 e T _ r b + _ R T b r r;i + _ (R b R i ) T r c;i = 0 e T v b + R T b ~ ! b r r;i + R T b R T i ~ ! i r c;i = 0 e T v b R T b ~ r r;i ! b R T b R T i ~ r c;i ! i = 0 e T v i : (3.17) where ~ is the skew symmetric matrix of coordinate vector, the entries of ~ are listed in Appendix A. The proof of the equality _ R T = R T ~ ! is available in [57]. With v b , ! b , v i and ! i dened, the derivative matrices are tabulated in Appendix A. For the MAV considered, the potential energy consists of the gravitational energy of each component and the elastic energy of each exure hinge, then we can have V = X i2b,1,2 m i gr T i u + X i21,2 1 2 k h;i 2 i ; (3.18) wherem i is the mass for theith component,u is the gravity direction vector,g is the 43 gravitational acceleration, k h;i is the stiness for the ith hinge. In the current case, u = [1; 0; 0] T andk h;1 =k h;2 . The matrix form of the generalized conservative force is @V @q T = X i2b,1,2 m i g @r i @q T u + X i21;2 k h;i @ i @q T i : (3.19) There are two types of non-conservative active forces for the system: the aerody- namic forces acting on the surfaces of the body and wings, and the torques driving the wings to ap and supporting the wings. These torques behave as active constraint forces resulting from the prescribed kinematic proles of the wings. Then the matrix form of the generalized non-conservative force is p q = @v b @ _ q T f ae,b + @! b @ _ q T ae,b 2 X i=1 cons;i ! + 2 X i=1 ( @v or;i @ _ q T f ae;i + @! i @ _ q T ae;i + R i cons;i ) (3.20) where f ae,b and ae,b are, respectively, the aerodynamic force and moment of the body evaluated at the center of the body. f ae;i and ae;i are the aerodynamic force and moment of the ith wing evaluated at the wing rooto r;i . The velocity v o r;i is the coordinate of the wing root velocityv o r;i in F 0 . We can nd from Eqn. (3.17) that v o r;i =v b R T b ~ r r;i ! b . The vector cons;i is the coordinate of the active constraint force cons;i in F b for the ith actuator. Substitute all the corresponding derivative matrices into Eqn. (3.7), @V @q T and (3.20), the equations of motion for the considered MAV is obtained as M(q) q + C(q; _ q)_ q +D(q; _ q) = 0; where D(q; _ q) = @V @q T p q ; (3.21) where M is the mass matrix, C is the structural damping matrix and D is the gener- alized external force vector. The rst 6 rows of Eqn. (3.21) represent the dynamics of the body, and the last 4 rows describe the motion of the wings. The dynamics of the body and wings are strongly coupled. In Eqn. (2.1), elements of q are assumed 44 independent of each other, but such independence is no longer valid in Eqn. (3.21) due to the introduction of constraints and constraint forces. These constraints must be handled directly or indirectly to make Eqn. (3.21) solvable. For a wing with fully prescribed kinematics, fully prescribed, the contributions of wings to body motion can be summed as an additional inertial term. However, when a wing is partially actuated with passive rotation mechanisms, the pitching motion must be solved together with the body motion, and the prescribed apping proles i (t) act as constraints to the MAV system. To account for this, the elimination method described in [39] was applied, which can satisfy the constraints automatically and also provide explicit formulas for the driving torques. In Eqn. (3.21), the unknown constraint forces cons;i are explicitly introduced. The apping motion rows must be removed from the original equations of motion to obtain a solvable reduced order system. This is the elimination method. The removal of the apping motion rows is straightforward for the inertial and potential force terms. The only unclear part is the non-conservative force term. We must understand the matrix structure of this term to separate the constraint forces from the other forces. Since the constraint forces only exist in the generalized force relating to the external moments, we dene the moment term p = @! b @ _ q T ae,b 2 X i=1 cons;i ! + 2 X i=1 @! i @ _ q T ae;i + R i cons;i =p ;1 +p ;2 ; (3.22) where p ;1 and p ;2 are the aerodynamic moment and constraint force components, respectively. p ;1 and p ;2 are p ;1 = 2 6 6 6 6 6 6 4 0 T T b ae,b + P 2 i=1 T T b R T i ae;i T T 1 ae;1 T T 2 ae;2 3 7 7 7 7 7 7 5 ; p ;2 = 2 6 6 6 6 6 6 4 0 0 T T 1 R 1 cons;1 T T 2 R 2 cons;2 3 7 7 7 7 7 7 5 : (3.23) 45 The part corresponding to the translational aerodynamic force is referred to as p f . The actuator drives the ith wing to ap around the x r;i axis and also sustains the bending induced by the wing. However, from the view of the entire system, the supporting torque to sustain the bending does no work when the deviation angle is always zero. The only active constraint force is the torque that drives the wing to ap, that is, cons;i = ;i e x r;i = b e T b [ ;i ; 0; 0] T = b e T b cons;i ; (3.24) where ;i is the magnitude of the driving torque for the ith wing, therefore we nd T T i R i cons;i = [ ;i ; 0] T ; p ;2 = [0 13 ; 0 13 ; ;1 ; 0; ;2 ; 0] T : (3.25) Substitute Eqn. (3.25) into Eqn. (3.21) and extract the apping motion rows, then we obtain the driving torques ;i = M (j;1:10) q + C (j;1:10) _ q + @V @q T (j) p f (j) p (j) : (3.26) where j = 7 + (i 1);i = 1; 2. The notation M (j;1:10) represents the submatrix of M with entries at the jth row and 1st-10th columns. The same rule applies to other matrices and coordinate vectors. When the apping rows are removed, the remaining rows in Eqn. (3.21) are the reduced order model without the constraint forces included M r (q) q r + C r (q; _ q)_ q r +D r (q) + 2 X i=1 A i i +B i _ i = 0; (3.27) where q r = [x b ;y b ;z b ; b ; b ; b ; 1 ; 2 ] T is the reduced coordinate vector. The matri- ces M r , C r , D r , A i and B i are available in Appendix A. Compared to Eqn. (3.21), Eqn. (3.27) can be easily transformed into a solvable state space representation. The elimination method can satisfy the constraints automatically, but we must 46 analyze the structure of the constraint force vector for dierent cases. For example, when the deviation angle of the wing is time-varying, the supporting torque becomes another active force. To handle dierent constraints in a unied way, the constraint- force method proposed in [57,60] is preferred. 3.4 Singularity-free and time-averaged formulations 3.4.1 Formulation with dual sets of Euler angles and quater- nions (a) (b) Fig. 3.2. The dual sets of Euler angles. (a) Coordinate frames describing the dual sets of Euler angles. (b) The relative rotation for the calculation of the roll angle 0 b . For free ight, the attitude of MAV can be in any physically allowed states. The pitch angle of the body can start from 0 and increase beyond=2. Using Euler angles in Sec. 3.1, the angular velocity of the body is restated as ! b = T b _ b = 2 6 6 6 4 1 0 sin b 0 cos b cos b sin b 0 sin b cos b cos b 3 7 7 7 5 2 6 6 6 4 _ _ _ 3 7 7 7 5 : (3.28) If b is =2, the matrix T b is singular. Since T b exists in the mass matrix M, the singularity can make the mass matrix singular. The fact behind the singularity ap- 47 pearance is that when the angular velocity is known, the yaw and roll motion of the body can be uniquely determined. Because the singularity only happens around b ==2, we employ two sets of Euler angles to resolve this issue. In Fig. 3.2(a), two inertial frames are dened to form the two sets of Euler angles. The inertial frame F 0 0 is obtained by rotating the frame F 0 about its y 0 axis for =2. The frame F b is the body-xed frame attached to the MAV body. The Euler angles b , b , b represent the rotation from F 0 to F b in the z-y-x convention. The auxiliary Euler angles 0 b , 0 b , 0 b represents the rotation fromF 0 0 toF b also in the z-y-x convention. Therefore, we can understand the two sets of Euler angles as consecutive rotations from dierent initial attitudes. FromF 0 , the body is initially vertical while, from F 0 0 , the body is initially horizontal. For both Euler angles sets, the nal step of the rotation is about thex b axis, thus the two inertial frames rotate to makex 0 or x 0 0 coincide with thex b axis, that is, 2 6 6 6 4 e x b e y 0 b e z 0 b 3 7 7 7 5 =R( b )R( b ) 2 6 6 6 4 e x 0 e y 0 e z 0 3 7 7 7 5 ; 2 6 6 6 4 e x b e y 00 b e z 00 b 3 7 7 7 5 =R( 0 b )R( 0 b ) 2 6 6 6 4 e x 0 0 e y 0 0 e z 0 0 3 7 7 7 5 : (3.29) We then obtain e x b = cos b cos b e x 0 + cos b sin b e y 0 sin b e z 0 (3.30a) e x b = cos 0 b cos 0 b e x 0 0 + cos 0 b sin 0 b e y 0 0 sin 0 b e z 0 0 = cos 0 b cos 0 b (e z 0 ) + cos 0 b sin 0 b e y 0 sin 0 b e x 0 (3.30b) The corresponding coordinates must be equal to each other, we can obtain the equal- ities below, 8 < : sin 0 b = cos b cos b cos 0 b = sin b =cos 0 b : (3.31) For the calculation of the roll angle b , the nal step of rotation is about thex b for 48 both Euler angle sets. Thus, the rotation axis is the same and we only compute the magnitude dierence of the rotation angle as in Fig. 3.2(b), b = 0 b b : (3.32) The axisy 0 b can be obtained through rotating axisy 00 b about the axisx b for , 2 6 6 6 4 e x b e y 0 b e z 0 b 3 7 7 7 5 =R() 2 6 6 6 4 e x b e y 00 b e z 00 b 3 7 7 7 5 : (3.33) With the coordinates for the axisx b , the equality is obtained sin b = sin b cos 0 b : (3.34) The switching algorithms are designed based on the Eqns. (3.31) and (3.34) for the two sets of Euler angles to avoid possible attitude singularities. Besides the dual sets of Euler angles method, the dynamic equations can also be formulated with quaternions. Because there are only two Euler angles for the description of the relative attitude of the wing with respect to the body, there is no singularity problem existing for the wing attitude representation. Here, we replace the rotational kinematic equation of Euler angles ! b = T b _ b (3.35) with the quaternion formulation [61] ^ ! b = 2 4 0 ! b 3 5 = 2E _ u; (3.36) 49 where u = [u 0 ;u 1 ;u 2 ;u 3 ] T and E = 2 6 4 u 0 u 1 u 2 u 3 u 1 u 0 u 3 u 2 u 2 u 3 u 0 u 1 u 3 u 2 u 1 u 0 3 7 5 : (3.37) The kinematics and corresponding derivatives are derived in the same manner as the Euler angles convention. The relations between the rotational matrices and quater- nion are dened in [61]. To resolve the unit constraintjuj = 1, the equations of motion formulated with Eqn. (3.21) by quaternions for body rotation and the relative Eu- ler angles for wing attitude are solved with the constraint-force method presented in [57, 62]. The elimination method is not available for the representation of Eqn. (3.21) with quaternion included. 3.4.2 Time-averaged models with mass-less wing assumption In addition to the singularity-free formulations, another important topic is the time- averaged model which is widely used in dynamics and controller design studies [9,63{ (c) (a) Stop Stop Body Motion Stop Stop Body Motion (b) Fig. 3.3. Computation work ows for various MAV dynamics models. 50 65]. If the rigid body assumption for the wing and the constant stiness model for the hinge are satised, the full dynamics model by Eqn. (3.27) can provide us with high- delity numerical results for the instantaneous dynamics of the MAV (Fig. 3.3(a)). However, due to the complexity of Eqn. (3.27), we cannot analytically solve them to obtain any closed-form solutions or carry out direct analyses to nd useful design criteria. For a MAV with fully prescribed wing kinematics, its full dynamics model can be reduced to a 6DOF rigid body model by neglecting the inertial contributions of the wings [27]. Here, we try to simplify our full dynamics model for the passive pitching case using the same approach (the small wing inertia assumption). M b = 2 4 I 33 m b 0 0 T T b J b T b 3 5 ; (3.38a) C b = 2 4 0 0 0 T T b J b _ T b _ T T b @T b @ b _ b 3 5 ; _ b = diag( _ b ;:::; _ b ) (3.38b) D b = 2 6 6 6 4 m b gu + X i2b;1;2 f ae;i T T b ae,b + X i21;2 T T b ~ r r;i R b f ae;i + T T b R T i ae;i 3 7 7 7 5 ; (3.38c) where I 33 is the identity matrix, ~ r r;i is the skew symmetric matrix of the wing root location vector (see Fig. 3.1(b)), M b , C b and D b are the corresponding matrices of M r , C r and D r with inertial contributions from the wings removed. Under the small wing inertia assumption, the coecients vectors A i and B i turn into zero directly. To solve the dynamic equations with the matrices in (3.38), the instantaneous pitching angles of the wings are still needed to compute the rotation matrix R i and the instantaneous aerodynamic forces f ae;i , ae;i . With the decoupled computation procedure shown in Fig. 3.3(b), for each time step, wing kinematics are solved rst using the body kinematics from the previous time step. Body motion is then solved using the new updated wing kinematics and Eqn. (3.38). To nd a simple and reliable method, we employ three strategies to estimate the wing pitching angle at each time 51 step. These are the full dynamics model (Model 1), the decoupled model (Model 2) and the single wing model (Model 3). In Model 1, the wing kinematics is updated using Eqn. (3.27). The updating of wing pitching angle is still coupled together with the calculation of body motion. In Model 2, the body kinematics is assumed to be xed when updating wing kinematics. The dynamic equation for the left wing is M (8;8) 1 + C (8;8) _ 1 = T T 1 ae,1 M (8;1:7) q r(1:7) C (8;1:7) _ q r(1:7) @V @q T (8) : (3.39) Note that ae;1 is the aerodynamic moment evaluated at the wing root, andq b =q r(1:6) and its derivatives are unchanged during a predictor-corrector numerical procedure. If we further neglect the inertial contributions from the body motion, the single wing model can be obtained as the Model 3. The equation of motion of the left wing is M (8;8) 1 + C (8;8) _ 1 = T T 1 ae,1 M (8;7) 1 C (8;7) _ 1 @V @q T (8) : (3.40) The 6DOF model of the body motion and the updating strategies of wing kinemat- ics are referred to together as the massless-wing models. Moreover, the time-averaged model can be obtained if the apping frequency of the wing is relatively much higher than the natural frequency of the body. The dynamic behavior of the body can be investigated at a time scale of one apping cycle, which can provide feasibility for long-endurance ight simulations. For the time-averaged model, the matrices M b and C b are in the same forms as those in Eqn. (3.38). The only change is the generalized external force term D b = 2 6 6 6 4 m b gu + X i2b;1;2 f ae;i T T b ae,b + X i21;2 T T b ~ r r;i R b f ae;i +T T b R T i ae;i 3 7 7 7 5 ; (3.41) in which is the operator for time-averaged value. Specially, the wing aerodynamic 52 moment is averaged in the body frame F b as R T i ae;i . As described in Fig. 3.3(c), the time step size of the inner loop is t in and that of the outer loop is t out . For each outer iteration to update the body motion, the inner loop runs for t out =t in times to compute the averaged aerodynamic forces and moments. During the updating process of the body motion, the wing kinematics is xed and the value is used for next outer loop iteration. The original purpose of developing massless-wing models is to simplify the full dynamics models so that we can carry out direct analyses. However, both the instan- taneous and averaged massless-wing models complicate the numerical simulation as shown in Fig. 3.3. The discussion of these massless-wing models in this approach is to validate whether we can make the small wing inertia simplication for the full dynamics models. If these models work, the next mission is to nd an alternative closed-form approximation for the wing kinematics instead of solving the instanta- neous dynamic equations. The results and discussions for the time-averaged model can be found in our publication [66]. 3.5 Integrated dynamics simulation with CFD solver and controller Under the weakly coupled FSI assumption, the ODE solver of Eqn. (3.27) is integrated with the Navier-Stokes equations solver to simulate the dynamic behavior of the entire MAV in free ight, as described in Fig. 3.4(b). In this block diagram, q pre represents the coordinates which are prescribed and q slv represents the coordinates to be solved from this integrated simulation. To be specic for the simulations in this section, the apping angles are assumed to be prescribed, i.e.,q pre = [ 1 (t); 2 (t)] T . The unknown coordinates are q slv =q r = [x b ;y b ;z b ; b ; b ; b ; 1 ; 2 ] T . On the ODE solver, to achieve a quasi-high order numerical scheme, the modi- ed Runge-Kutta scheme with extrapolated aerodynamic forces are employed. Eqn. 53 (a) (d) (b) A A C B D (c) Fig. 3.4. Integrated dynamics simulations with CFD solver and controller. (a) Integrated simulation for the entire MAV. (b) Mesh topology for the CFD simulation of the entire MAV. (c) One slice for the assembled grids. (d) Simulation with controller included. (3.27) is rst re-formulated as 8 < : _ q r = _ = q r ;;q pre ; _ q pre ; q pre ;f ae ; ae : (3.42) During the numerical solving procedure, the aerodynamic forces f ae and torques ae from CFD solver are discrete values without explicit formulas. The r-order poly- nomial functions ^ f ae (r;N;t) and ^ ae (r;N;t) are constructed with least-square tting from the latestN points extrapolate for the prediction of aerodynamic forces in next time steps. The equation is re-expressed as ^ _ = q r ;;q pre (t); _ q pre (t); q pre (t); ^ f ae (r;N;t); ^ ae (r;N;t) : (3.43) 54 The standard 4th order Runge-Kutta procedure were applied to Eqn. (3.42) in our simulations. These above steps constitute our modied Runge-Kutta scheme with extrapolation. Since the extrapolation operations can introduce numerical errors, the time step size should be small enough to guarantee the accuracy of extrapolation. Due to the essence of prediction-correction process, the modied Runge-Kutta scheme is expected to be in higher order than the Euler scheme and permits larger time step size than the Euler scheme, which was observed in our simulations. On the Navier-Stokes solver, the aerodynamic forces are computed in commercial package STAR-CCM+ using the overset grid approach. The mesh topology for the CFD simulations is depicted in Fig. 3.4(b). L1 is the mesh block for the background air ow. L3 and L4 are the mesh blocks for the body. L5 and L6 are the mesh blocks for the two wings, respectively. L2 is the transition block as a container of the two body mesh blocks. The interface between L1 and L2 is in overset grid type. The interface between L2 and L3 (or L4) is in interior type. The interface between L5 and L3 is also in overset grid type (the same for L6 and L4). The interface between L3 and L4 is in interior type. During the CFD simulation, the background mesh is xed in space. The block L2 moves with the translational and rotational velocities of the body (center of mass). The body meshes are relatively xed to L2. L5 moves with respect to L3, and L6 moves with respect to L4. To guarantee the successful assembling of the overset grids, there is a gap between the wing root and the body surface. The gap can be treated as the length of the end eector part which is outside the body surface. The mesh domain for the body is split into two blocks (L 3 and L 4 ) to ensure symmetric mesh assembly. One slice of the assembled 3D grids depicting the above topology is visualized in Fig. 3.4(c). High-resolution grids are provided for the wing and the boundary layer region of the body. Alternative topologies for the wing mesh are discussed in Sec. 2.2.2. Because the mesh block for the background is too large, only part is drawn in this gure. Inside STAR-CCM+, several steps must be completed to guarantee the 55 correct conguration of the mesh assembly. A straightforward application of the dynamics model is to treat it as a physical representation of the MAV in the closed-loop simulation with a controller included, as shown in Fig. 3.4(d). Here, the control diagram is adapted from [67] with the robotic plant replaced with the full dynamics model (block A). Since, in our current study, the actuation signals are the apping angles instead of the electrical signals, the output from the controller block (C) is u p = [ 1 ; 2 ] T . The actuation block (B) is a constant value 1. The signal of the apping angle is in a sinusoidal form as i (t) =a i sin(2t) +b i (t); (3.44) If the apping frequency is xed during ight, a i (t) and b i (t) are the true control signals. Inspired by the control algorithms in [67], a simpler PID algorithm for the updating of a i (t) and b i (t) are employed as a i (t) = 0 (t) + K a,p e b + K a,i Z e b dt + K a,d _ e b ; (3.45a) b i (t) = 0 (t) + K b,p e b + K b,i Z e b dt + K b,d _ e b (3.45b) where 0 is the baseline amplitude, 0 is the initial bias,e b =q b q b,ref is the relative error of the body states with respect to the reference, and K a,p , K a,i , K a,d , K b,p , K b,i and K b,d are the PID parameters for our controller. In our simulations, non-zero K a,p , K a,i and K a,d are for the height and roll control of the robot, and non-zero K b,p , K b,i and K b,d are for the horizontal drifting and pitch motion of the body. The mechanism of force generation for these schemes are discussed in the controller design of [67,68]. For the control law in 3.45a and 3.45b, the values of i are not enough for our dynamics simulation. The time rate and acceleration of apping angles must be estimated. In our simulation implementation, the history values of a i and b i are saved. Here, the control signals are updated in a larger sampling time T than 56 t > t n STAR-CCM+ D 1 D 2 D m User Code (C++) STAR-CCM+ D 1 D 2 D m User Code (C++) t t MATLAB MATLAB No FILE MANAGER output.dat q n read File locking Dyn Cont Dyn Cont Yes write (loop until t = t n ) q n q n Fig. 3.5. Parallel implementation of the integrated simulation with STAR-CCM+ and MATLAB. numerical step t). The spline-based approximations of ^ a i (r;N;t) and ^ b i (r;N;t) are constructed to nd approximations for _ a i , a i , _ b i and b i . Therefore, the angular velocity and acceleration of apping angle are, respectively, _ i (t) = _ ^ a i sin(2t) + ^ a i (t)(2) cos(2t) + _ ^ b i ; (3.46a) i (t) = ^ a i sin(2t) + 2 _ ^ a i (t)(2) cos(2t) ^ a i (2) 2 sin(2t) + ^ b i : (3.46b) The variations ofa i andb i are in low frequency comparing to the apping motion. The terms ^ a i (t)(2f) cos(2ft) and^ a i (2f) 2 sin(2ft) dominate the values of _ i and i , respectively. In our study, the purpose of closed-loop simulations is to exemplify the eectiveness of this framework. For the control of real apping-wing MAVs, the algorithms in [18, 67, 68] can be implemented. Disturbance models can also be introduced as in block B. The coding of the integrated simulation for a single wing in serial model is simple, only the C++ \User Code" interface of STAR-CCM+ is required to integrate the 57 a (a) (b) Fig. 3.6. Blade element analyses for the wing and the body in free ight. (a) A wing with an arbitrary wing root kinematics. (b) Body motion and wing induced air ow. dynamics equation (2.7) to the simulation. Unfortunately, the \User Code" inter- face does not provide a straightforward mechanism for the parallel implementation, especially for the case in which les are employed for data exchange and communi- cation. Fig. 3.5 describes the implementation of the parallel interface code for the full dynamics simulation in which the dynamics equations (3.27) and controller code are solved or updated with MATLAB. The C++ \User Code" and MATLAB com- municate loosely with the data le `output.dat' for the generated velocity vector q n . To realize the writing of `output.dat' correctly in the parallel mode, the le locking mechanism (Windows or Linux le system) are employed. The code for le locking and call of MATLAB are listed in the Appendix C. 3.6 Quasi-steady models for a MAV in free ight 3.6.1 Wing with arbitrary root kinematics To estimate aerodynamic forces acting on a wing during free ight, wing root kinemat- ics must be considered when computing the characteristic velocity and acceleration for each wing section. Since the existence of wing root kinematics is possibly able to change the spatial structure of the ow eld around the wing [21,69], the parameters calibrated for hovering ight might be inaccurate for free ight. To extend the models for hovering ight to free ight, special conditions or corrections must be satised or 58 specied, respectively. As depicted in Fig. 3.6(a), stroke motion frame F s;i fo r;i ,x s;i ,y s;i ,z s;i g is intro- duced to dene characteristic angles and kinematics for each wing section. To be specic,x s;i coincides withx r;i , andy s;i is the carried axis ofy r;i after a rotation of i about x s;i . For each wing section, characteristic velocity is the incident velocity v n inside thex s;i z s;i plane at the intersection point between the spanwise axis of the hinge and the wing section itself. v n is computed as v n = v n;z s;i e z s;i +v n;x s;i e x s;i ; (3.47a) v n;z s;i = v r e z s;i + ! i e x s;i r; v n;x s;i = v r e x s;i ; (3.47b) v r = b e T R b v b + ~ ! b r r;i : (3.47c) wherev n;z s;i andv n;x s;i are the components alonge z s;i ande x s;i axes, respectively, and v r is the wing root velocity. v r includes contributions from both the body translational motion v b and the body rotation ! b . The angle of attack is computed as (r) = hv n ;e x i i with the operatorh;i dening the angle between two vectors. For the calculation of f rot,1 , replace _ with s;i ! i(2) in Eqn. (2.13b). Since f rot,2 represents the damping eects from the rotation of wing relative to the body, the formula is the same as (2.13c) for a single wing. To estimate the added mass force, the characteristic acceleration is still the one at the midpoint of each wing section. The normal component of the characteristic acceleration is a a,n =a a e z i ; (3.48a) a a =a r + r a =a r + 0 e T (R i R b ) T ~ ! i ~ ! i + _ ~ ! i r a ; (3.48b) a r = 0 e T _ v b + _ R T b ~ ! b r r;i + R T b _ ~ ! b r r;i : (3.48c) wherer a is the position vector from wing root to the midpoint of wing section under consideration. The accelerationa r is the time derivative ofv r . Direct superposition in (3.47) and (3.48) is consistent with the handling of wing 59 root kinematics in [28, 34, 66]. As discussed in [21, 69, 70], wing velocity is to change the ow structure around the wing compared to the hovering case. Accuracy of the direct superposition can be aected by magnitudes of wing root velocity and acceleration. The corresponding aerodynamic coecients for the hovering case might also not be applicable to the free ight cases. However, there lacks generalized and reliable models for complex wing root kinematics. To make the model more applicable and simple, lumped correction coecients are introduced for the cases with large wing root velocity and accelerations as ^ v r = v v r ; ^ a r = a a r ; (3.49) where v and a are the correction coecients for velocity and acceleration, respec- tively. v and a account for the nonlinear superposition of the wing root kinematics. v primarily re ects the corrections resulting from two aspects: (1) The ow struc- ture around the wing is changed by the incoming ow; (2) Pitching prole of the wing is possibly altered because of body rotations. a re ects the violation of the simple formula for the equivalent volume of added mass in Fig. 2.4(d). For a condi- tionally superposition model, ^ v r and ^ a r are to replacev r anda r in Eqns. (3.47) and (3.48), respectively. Conditions for the direct superposition and possible corrections are discussed in Sec. 3.8.3, Sec. 3.8.4 and Sec. 3.8.6. 3.6.2 Simplied cylindrical body For a MAV in free ight, wings dominate aerodynamic force generation on the sys- tem [28]. However, without full dynamics simulation, it is unknown how signicantly the aerodynamic force on the body can change the motion of the MAV. To quan- titatively investigate the accuracy of the derived full dynamics model, the inclusion of aerodynamic force models for the body is crucial. When moving in space, body also lies inside the air ow induced by the wing. If the apping frequency of the wing 60 is high enough, the induced velocity can be comparable to the translational velocity of the body. As seen in Fig. 3.6(b), the variation of the body cross section along the cylinder axis is assumed to be small, so we only consider the force component normal to the cylinder axis. For each element dl , the Morison equation is employed to estimate the drag force on each element [71] df b,d = 1 2 C b,d jv f,n jv f,n d e dl + 1 4 C b,a a f,n d e 2 dl; (3.50a) v f,n =v b,n +! b l e +v ind ; (3.50b) a f,n =a b,n + b e T ~ ! b ~ ! b + _ ~ ! b l e + _ v ind ; (3.50c) On the right-hand side of Eqn. (3.50a), the rst term is the translational force and the second term is the added mass force. C b,d is the drag coecient and C b,a is the added mass correction coecient. v f,n is the normal component of the body velocity relative to the air anda f,n is the normal acceleration. Here, a characteristic diameter d e is used instead of the local value to simplify calculation. In Eqn. (3.50b), v b,n is thez b component of the body translational velocity,l e is the location vector of the element andv ind is thez b component of the wing induced velocity. On acceleration, a b,n = b e T a b(3) is the z b component of the body translational acceleration. The wings rotate together with the body, so the time rate of the velocity magnitude v ind is not necessarily the true acceleration magnitude of the induced air ow. To simplify calculation, the time rate of v ind is directly treated as the magnitude of the acceleration, i.e. _ v ind = b e T _ v ind . To approach practical calculation, we must rst nd a way to estimate the induced air ow velocity. Rather than examining the details of the ow eld, we propose the existence of a critical point on the wing along the hinge axis atr =r c . The velocity at this point is the characteristic velocity of the entire induced air. Thez b component of the critical velocity and the induced velocity distribution are, respectively, v ind,c = _ mod r c cos c cos; mod = (t t ); (3.51a) 61 v ind =v ind,c 1c loc l r l b ; (3.51b) where c is the critical pitching angle of the wing corresponding to the maximum apping angular velocity, mod is the modied apping angle with a phase shift and the corresponding shift time is t . The induced air ow still satises the no-slip condition and the cos term represents the velocity component normal to the wing. v ind,c is the value at the critical point around the wing root, and a linear pattern is assumed to model the local induced velocity along the body axis, where l b is the length of the body, l r is the distance from the element to the wing root, and c loc is the slope. The parameterc loc is adjusted according to experimental or CFD data for a specic body shape. In our study, we only addresses the longitudinal motion. The direction of v ind is always along the z b axis. For a MAV with two wings, a direct superposition of the induced velocity is assumed as v ind, t =v ind,1 +v ind,2 : (3.52) For the simulations of a robotic MAV immersed in water in Sec. 3.6, the wing induced ow is neglected for modeling simplicity. For a robotic MAV with complex airframe (e.g., the one in Fig. 3.7(a)), a characteristic diameter d e is rst picked for the body based on the geometry and CFD simulations are then employed to nd the QS parameters, such as C b,d , C b,a . Besides drag due to motion, static buoyancy force on the body is also considered. Unsteady pressure contributions such as the Froude-Krylov force is not modeled for simplicity. The vertical total buoyancy force is estimated as f buoy =V disp g; (3.53) where V disp is the volume displaced by the body. For the case in light media (e. g., air), the buoyancy force on the body is expected to be very small. 62 (a) (c) (b) Fig. 3.7. The 70-mg insect-sized robotic MAV and related experimental settings. (a) CAD model for the MAV II. (b) Front view of the open-loop test in the air. (c) Front view of the open-loop test in water tank. To implement the integrated simulation with the QS models, the CFD solver in 3.4(a) is replaced with the QS models for the wing and body. The QS parameters for the wing are calibrated with the CFD simulation for a single wing. More advanced calibration method is discussed in 3.8.3. If the induced air ow is not considered, the QS parameters for the body are obtained with the incoming ow simulation or sinusoidal oscillation. 3.7 Preliminary experimental validation 3.7.1 Options for experimental validation To validate the proposed full dynamics model, a 70-mg robotic MAV (MAV II) is de- signed as shown in Fig. 3.7(a) and fabricated with the method explained in the next subsection. The supplemental video for the experimental tests and associated simula- tions is available at https://www.uscamsl.com/resources/dissertationchang/S1.mp4 and https://www.metasimulating.com/resources/dissertationchang/S1.mp4. In the later sec- tions, the supplemental video is referred to as S1.mp4. For the purposes of quantita- tive or qualitative validation, three options are available for experimental tests: 1. Open-loop jump test in the air in Fig. 3.7(b). In the jump tests, the apping proles of the two wings are nominally symmetric to each other about the x b z b plane. The MAV typically jumps from the ground and quickly crashes back (see the supplemental S1.mp4). The ight time in the air is very short for most of our tests. 63 Even for a long enough ight, the trajectory of the MAV also quickly changes from simple 2D to an unexpected 3D curve. During the tests, the MAV is not perfectly fabricated. The legs have strong interactions with the ground before taking o so that the initial conditions are aected. There might also exist air or wire disturbances. Since the MAV is inherently unstable [66], it is challenging to achieve repeatable and quasi-2D tests for quantitative validation. Most importantly, for the two-winged conguration in Fig. 3.7(a), the yaw motion is not eectively damped [16,18,72,73], small fabrication errors and asymmetric initial conditions can trigger strong yaw rotations. With the above observations and current fabrication quality for the MAV, it is not practical to quantitatively validate our proposed full dynamics model with open-loop tests in the air. 2. Closed-loop ight test in the air. Ideally, repeatable results could be expected to acquire in fully controlled experiments, in which disturbances from fabrication errors and case-varying initial conditions are completely rejected. However, the design conguration of the MAV in Fig. 3.7(a) is not eective in damping yaw rotations and resisting horizontal drifting, as shown in the supplemental S1.mp4. Even with a controller included [18], the tests are not quantitatively repeatable. The above reasoning is consistent with the results observed in [34], in which bad quantitative matching is seen for a control test. From another view, in the closed-loop numerical simulation, the controller can compensate for modeling errors. The accuracy of the full dynamics model cannot be addressed when a controller is presented. Thus, the closed-loop test in the air is also not preferred for quantitative validation. 3. Open-loop jump test in the water tank. Inspired by the results in [6] and [74], the apping-wing system is expected to be more damped in a heavier media than air. After implementing the underwater tests as in Fig. 3.7(c), quantitatively repeatable and quasi-2D tests were observed for the apping frequency at 13 Hz. With the quasi- 2D tests, there are fewer concerns on how to exact the real trajectory of the MAV, but more challenges come from extra modeling corrections. The required modications 64 are described in Sec. 3.7.3. After the above analysis, the underwater tests are preferred for the quantitative validation of our full dynamics model. Qualitative validation with open-loop tests in the air is also presented to demonstrate the capability of the model in capturing the in-air ight physics of the considered MAV. 3.7.2 Experimental setups for underwater and in-air tests The prototype in this paper was fabricated by employing the smart composite mi- crostructure (SCM) method in [17,18,75]. The overall process mainly includes three steps. First, bond multiple layers of composites with dierent mechanical properties, e.g. carbon ber/epoxy composite and polyester lm for the wing, into a laminate using adhesive lms (Pyralux FR1500, DuPont); then precisely cut individual compo- nents from the laminates by utilizing a diode-pumped solid-state (DPSS) laser (DCH- 355-3, Photonics Industries) with a spot diameter of 10m; lastly manually assembly all components, including pieces of airframe, transmissions, exural hinges and wings with their mating features under a optical microscope, and apply cyanoacrylate (CA) glue (Loctite 416, Henkel) to strengthen their adhesion. We also exploited the pre- stack method introduced in [76] to fabricate the PZT bimorph actuators for the MAV. The completed prototype is shown in Fig. 3.7(a). To enable operating in an aquatic environment, the robotic MAV was electri- cally isolated by applying a thin layer of CA glue on the surfaces and bases of PZT actuators, where connections of tethered wires locate. To avoid any possible electri- cal shorting, deionized water (Resistivity 18 megohms/cm, Ecoxall Chemicals) was chosen as the ow media. The driving voltage240V is applied to the actuators. 3.7.3 Modeling modications for the underwater tests To achieve quantitative validation with repeatable and quasi-2D experimental results, underwater tests were carried out for our robotic MAV as the setup in Fig. 3.7(c) 65 and the clip C3 in the supplemental S1.mp4. However, the full dynamics model in Sec. 2 derived for a perfect MAV in the air cannot be directly applied to a robotic MAV moving in the water. Specic model modications or kinematics adjustments are made to produce reasonable simulation results for our underwater tests. In Eqns. (3.23) and (3.27) for a MAV in the air, aerodynamic forces and torques are lumped quantities. The numerical procedure for (3.27) denes a numerical weakly coupled FSI procedure as in [23] and [66]. However, for the underwater tests, the weakly coupling assumption is not valid anymore because the added mass force is too large. If the added mass force is estimated in the lumped form on the right-hand side of the discretized equation, the one-time-step delay can scale up the numerical errors into innity in limited steps. Converged results cannot be achieved with the direct discretization of Eqn. (3.27). Thus, the added mass terms corresponding to q must be extracted and added into the two sides of Eqn. (3.21). A more strongly coupled form is formulated as (M + M a ) q =C_ qDp a ; (3.54) p a =M a q; (3.55) where p a is part of the total added mass force, which is directly associated with q. M a is the added mass inertia matrix. The corresponding reduced form can be derived similarly to that of Eqn. (3.21). For the full dynamics model of a MAV in the air, the apping prole is prescribed. To drive the wings to ap, sinusoidal electrical signal with constant bias is sent to the PZT actuator as e(t) =e 0 sin(2t) +e ; (3.56) where e 0 is the amplitude of the electric signal, is the signal frequency and e is a constant signal drifting. The mapping between the apping angle and electrical 66 signal can be described in a linear and time-invariant model [38]. Accordingly, for asymmetric electrical signal, the apping prole also has a constant bias. However, as observed in the supplemental S1.mp4, the constant bias model for the apping prole is not valid for the underwater tests. Strong interactions among water, wing, transmission and the actuator can change the dynamic response of the actuator. For simplicity, a piecewise time-varying drifting model for the apping prole is assumed to exist for the tests in Sec. 4.1 as (t) = 0 sin(2t) + 0 (t); (3.57) 0 (t) = 8 < : a 0 t +b 0 tt c a 0 t c +b 0 t>t c ; (3.58) where 0 is the amplitude of the apping prole, a 0 is the drifting speed, b 0 is the initial value of the drifting. 0 becomes constant after t c . As the mapping between the electric and apping signals is not modeled, the values of parameters for are estimated from the test video. In Sec. 2.2.1 and Sec. 3.5, the exure hinge is assumed to pitch freely. However, there possibly exists geometric collisions between the upper and bottom shoulders [23], which places a geometric limit to the pitching angle. After the pitching angle is limited by hinge geometry, the wing can deform to achieve a larger pitch angle. Modeling this kind of elastic one-side constraint is not easy. For simplicity, we employ the concept of extra bending torque to model the comprehensive eect by hinge geometry limit and wing deformation, extra;i =k c ( i limit ) +c c _ ; i > 0 > 0 (3.59) wherek c is the correction stiness,c c is the correction damping, and limit is the angle limit. A negative sign is applied for the case i < 0. The extra torque is added into the row corresponding to the elastic torque of hinge i in the equations of motion 67 (3.27). The values of k c and c c are estimated based on the amplitude of the pitching angle of the wings in the clip of supplemental video S1.mp4 for Sec. 4.1. During the initializing of a jump test, the legs of the MAV interact with the ground before taking o. Extra elastic bouncing force and torque are estimated for the contacting point between the leg and the ground for x leg <x leg (0) as f b =k leg (x leg x leg (0)); (3.60a) b =k leg (x leg x leg (0))z leg ; (3.60b) where k leg is the stiness of the leg, (x leg , z leg ) denes the location of the contacting point of the front/rear leg with the ground, andx leg (0) is the initial vertical position of the leg tip att = 0. To void more complicated data processing and model adjustments, the eects of the water refraction, the power wire and the interaction between water wave and tank wall are not modeled. Based on our observations during the tests, these eects are not important for the quasi-2D results. 3.8 Results 3.8.1 Integrated simulation with CFD solver and controller The integrated simulations described in Fig. 3.4 were implemented for the open- loop and closed-loop ight of MAV I. The geometric and inertial parameters for the body, wing and hinge of MAV I are listed in the Appendix B. For the open-loop simulation, the apping prole is Eqn. (3.44) with a i (t) = 0 and b i (t) = 0 = 0. The states of the MAV are initially disturbed by the non-zero apping velocities. For the closed-loop simulation, the input signal is also Eqn. (3.44) with the PID prole (3.45). The controlled state is b and the reference value is b = 0. For the controller, only the proportional block (P block) is employed, K b,p(5) =5. The proles of the input signals are recorded in Fig. 3.8(a). As seen in Fig. 3.8(b)-(c), compared to 68 Pitch angle of the body 0 0.4 0.8 1.2 Time (s) 0 0.04 0.08 0.12 0.16 0.2 θ b (rad) closed-loop x b (m) z b (m) 0.04 0.08 0.12 0 0.02 0.04 0.06 closed-loop Planar position and attitude 0 open-loop Flapping angles for open-loop and closed-loop simulations Time (s) 0 0.04 0.08 0.12 0.16 0.2 -2 -1 0 1 2 closed-loop (rad) (a) (b) (c) Fig. 3.8. Integrated dynamics simulations with CFD solver for MAV I. (a) Instanta- neous apping angles. (b) Planar position and attitude. (c) Body pitch angle. Note: the open-loop simulation terminated around 0.19 s because the moving mesh reaches the boundary of the background mesh. 69 (a) (b) (c) (d) Fig. 3.9. Flow elds for the open-loop ight of the MAV. (a)-(b) Pressure and velocity magnitude elds for the longitudinal plane of the body. (c)-(d) Pressure and velocity magnitude elds for the wing section with maximum chord length. the open-loop ight, the pitch angle of the body with a controller is conned in a smaller range and the horizontal drifting also becomes signicantly smaller. Though simple, the controller is very eective for this MAV system. If more disturbances exist, the structure of the controller can be more complicated, but it does not aect the structure and implementation of the integrated simulation. For the open-loop simulation, the pressure and velocity elds around the body are visualized for the region within the L2 layer (see Fig. 3.4(c)) in Fig. 3.9(a)-(b). It is seen that the mesh for L1 is very coarse, which can signicantly reduce the computing burden. In Fig. 3.9(c)-(d), the pressure and velocity magnitude elds are visualized for the region within the L3 layer. If comparing the visualized results, the ow eld behind the body is more complicated than that around the wing. The magnitude of the pressure behind the wing is around ten times larger than that behind the body. From an intuitive view, the aerodynamic forces on the wing should be much larger 70 than these on the body. The magnitudes of the velocities are consistent with the results of the pressure distribution. The visualization of the ow eld can be useful in analysis when dierent wing shapes or kinematics are applied for the MAV. The integrated simulations with CFD depend on very few empirical parameters, which is an ideal property for the parametric design of the apping-wing MAV with dierent geometry or kinematics. However, problems appear when computing e- ciency is considered. For the simulation of a 0:2 s physical ight, the simulations take around 3 4 days on 64 CPU cores on the USC HPCC cluster. If more CPU cores are employed, the computing eciency becomes even lower. If the MAV ies in a larger space or a longer physical time, the required CPU hours can be tremendously huge. The integrated simulation with STAR-CCM+ does not allow real massive par- allel computing. Since the exact implementation of the overset grid technique in STAR-CCM+ is unknown, the bottleneck of computing eciency prevents the inte- grated simulations from acting as a useful tool for fast iterative mechanical design or controller tuning. Also, the grid assembly failure happens frequently for our mesh settings. Our proposed fast overset grid assembly algorithm is discussed in Chapter 4 for our in-house CFD solver. 3.8.2 Adaptability of the quasi-steady models for the wing and body in free ight The adaptability of the QS models for a wing in free ight was tested with dierent incoming ow velocities. In Fig. 3.10, simulations were performed for U = 0:5 m/s (Figs. (a)-(c)) and U = 1:0 m/s (Figs. (d)-(f)) with parameter set P 2 in Table 2.3. Since the ow eld is initialized with zero velocity, the transition procedure can be observed in the rst few apping cycles. For the U = 0:5 m/s case, there is a satisfactory agreement between the quasi-steady and CFD results. However, for the U = 1:0 m/s case, large discrepancies exist in the lift curves around the peak 71 CFD QS-P 2 Time (s) 0 0.02 0.04 0.06 0.08 0.1 0 1 2 -1 0 1 -1 0 1 0 1 2 -1 0 1 -1 0 1 Lift for Δθ t =0.5π and U = 0.5 m/s (a) (b) (c) (d) (e) (f) Time (s) 0 0.02 0.04 0.06 0.08 0.1 Time (s) 0 0.02 0.04 0.06 0.08 0.1 Time (s) 0 0.02 0.04 0.06 0.08 0.1 Lift for Δθ t =0.5π and U = 1.0 m/s Drag for Δθ t =0.5π and U = 0.5 m/s Drag for Δθ t =0.5π and U = 1.0 m/s Pitching moment for Δ θ t =0.5 π and U = 0.5 m/s Pitching moment for Δθ t =0.5π and U = 1.0 m/s Time (s) 0 0.02 0.04 0.06 0.08 0.1 Time (s) 0 0.02 0.04 0.06 0.08 0.1 Fig. 3.10. CFD validation of the quasi-steady model for the case with incoming air ow. (a)-(c) Lift, drag, pitching moment for U = 0:5m/s. (d)-(f) Lift, drag, pitching moment for U = 1:0m/s. or valley regions. Since, for both cases, the incident velocity is calculated with the same formula (3.47), the discrepancies in Figs. 3.10(d)-(f) imply that the validity of the blade element analysis can deteriorate for a large velocity. Also, the CFD results present stronger variations on the amplitudes of the lift curves for larger incoming ow velocity. The incoming ow can strengthen the unsteadiness of the ow. If the forward ight velocity increases further over 1:0 m/s, it is possible that the quasi- steady assumption cannot be satised for the current model. Accordingly, the direct superposition formula Eqn. (3.47b) should be corrected for large incoming velocities, e.g., using the lumped correction Eqn. (3.49). 72 CFD validation of the quasi-steady model for the cylindrical body is plotted in Fig. 3.11. Since the Morison equation has been widely investigated for the ow around cylindrical shape [71], in this subsection, we only validate the proposed model for the body force resulting from wing induced air ow. Comprehensive simulations with the Morison equation included are presented for the entire MAV in Sec. 3.8.4. The body is xed in space and the two wings ap with respect to the body. The values of 0 and 0 are the same as those for Fig. 3.10. The corresponding parameter sets are listed in Table 3.1. The critical pitching angle is calculated with c = 0 sin( t ). Eqn. (3.51) represents an averaging eect of the air ow around the body, so it is challenging to nd the specic local velocity to compare. Thus the reliability of the model is indirectly demonstrated by the accuracy of aerodynamic force estimation. Time (s) 0 0.01 0.02 0.03 0.04 0.05 -0.05 0 0.05 CFD QS Time (s) 0 0.01 0.02 0.03 0.04 0.05 -0.5 0 0.5 Time (s) 0 0.01 0.02 0.03 0.04 0.05 -0.1 0 0.1 Time (s) 0 0.02 0.04 0.06 0.08 0.1 -0.05 0 0.05 Time (s) 0 0.02 0.04 0.06 0.08 0.1 -0.2 0 0.2 (a) Time (s) 0 0.01 0.02 0.03 0.04 0.05 -0.2 0 0.2 (b) (c) (d) (e) (f) Fig. 3.11. CFD validation of the quasi-steady model for the body force resulting from the wing induced air ow. 73 Table 3.1: Parameters for the quasi-steady model of the body. Name t C b,d C b,a c loc r c (mm) 2t c d e (mm) P 6 0:5 50 Hz 1.0 1.0 1.8 1.13 0:25 0.995 3 P 7 0:5 100 Hz 1.0 1.0 1.8 1.25 0:25 0.995 3 P 8 50 Hz 1.0 1.0 1.8 1.13 0:00 0.000 3 P 9 100 Hz 1.0 1.0 1.8 1.25 0:00 0.000 3 P 10 0:75 100 Hz 1.0 1.0 1.8 1.25 0:15 0.704 3 Results for the single frequency cases (100 Hz) are shown in Figs. 3.11(a)-(b), (d)- (e). For both the t = 0:5 and t = cases, the proposed quasi-steady model is able to capture the basic characteristics of the aerodynamic force generation on the body. Fig. 3.11(c) and (f) show the mixed frequency case for t =. The left wing aps at 100 Hz and the right wing aps at 50 Hz. The results indicate that the proposed model performs consistently well for the mixed frequency case, which highlights the applicability of the direct superposition model (3.52). The above results also re ect the accuracy of the Morison equations. However, as shown in Fig. 3.11, the t = case matches better than the t = 0:5 case, indicating the accuracy of the proposed model is not consistent for dierent wing kinematics. Considering the complexity of the induced air ow, the estimation accuracy in Fig. 3.11 is adequate as a rst approximation to study the role of body aerodynamic forces in MAV ight. The calibrated parameters in Table 3.1 are used in the next section. 3.8.3 Full dynamics simulation with quasi-steady models The integrated simulations with CFD solver and QS models share the same inertial and conservative forces parts in the equations of motion (3.27). The full dynamics simulation with a CFD solver can directly verify the accuracy of the QS models for free ight. Integrated simulations were implemented for MAV I with two kinematic condi- 74 x b , z b (m) Time (s) 0 0.04 0.08 0.12 0.16 0.2 -0.15 -0.10 -0.05 0 0.05 θ b (rad) Time (s) Pitch angle of the body 0 0.04 0.08 0.12 0.16 0.2 -0.5 0 0.5 1.0 Time (s) Pitch angle of the body 0 0.04 0.08 0.12 0.16 0.2 -2 -1 0 1 Planar coordinates of the body centroid θ b (rad) x b -QS-P4 (Crot = 1.0) x b -CFD ( 0 = 0.0) x b -QS-P5 (Crot = 1.6) z b -QS-P4 (Crot = 1.0) z b -CFD ( 0 = 0.0) z b -QS-P5 (Crot = 1.6) QS-P4 (Crot = 1.0) CFD ( 0 = 0.0) QS-P5 (Crot = 1.6) QS ( 0 = 0.1, P2) CFD ( 0 = 0.1) QS ( 0 = 0.0, P2) CFD ( 0 = 0.0) (a) (b) (c) Fig. 3.12. Comparisons of the body states between the CFD and QS results for MAV I. (a) Planar coordinates of the body centroid for 0 = 0. (b) Body pitch angle for 0 = 0. (c) Body pitch angle for 0 = 0:1. 75 f 1,1 (mN) Time (s) Vertical force on the wing 0 0.04 0.08 0.12 0.16 0.2 -1 0 1 2 (a) f 1,3 (mN) Time (s) Horizontal force on the wing 0 0.04 0.08 0.12 0.16 0.2 -2 0 2 CFD ( 0 = 0.0) QS-P5 (Crot = 1.6) (b) τ 1,2 (µNm) Time (s) Pitching moment on the wing about its hinge axis 0 0.04 0.08 0.12 0.16 0.2 -1 0 1 (c) Fig. 3.13. Comparisons of the instantaneous forces on the left wing between the CFD and QS results. (a) Instantaneous vertical force in the inertial frame. (b) Instantaneous horizontal force in the inertial frame. (c) Pitching moment of the wing about the hinge axis in wing-xed frame. 76 tions: (I) 0 ==3, 0 = 0:0 and (II) 0 ==3, 0 = 0:1 with both CFD solver and QS models. The simulations with CFD solver share the same settings with the results in Sec. 3.8.1. As seen in Fig. 3.12(a)(b), the dynamic behavior of the MAV generated from the QS models resembles that of the CFD solver, which directly demonstrates the eectiveness of the QS models and indicates that both sets of parameters re ect the true physics of the MAV dynamics. However, in terms of accuracy, the pre- cisely calibrated parameter set P4 provides bad quantitative matching with the CFD results. 50% relative errors can be seen for the maximum body pitch angle. The superposition correction c was tried, but it did no help on the prole of the body pitch angle. Since b can directly aect the rotational components in instantaneous force generation, parameter set P5 with larger C rot signicantly improved the accu- racy of the QS model for all body states seen in Fig. 3.12(a), (b). To further show the prediction of the parameter set P5, the simulated pitch angle of the body with 0 = 0:1 is plotted in Fig. 3.12(c). For both 0 s, consistent accuracy is achieved for the parameter set P5, which indicates C rot = 1:6 re ects the contributions from the rst rotational component than the value of 1.0. The instantaneous vertical and horizontal wing forces in the inertial frame and the pitching moment of the wing in the wing-xed frame are plotted in Fig. 3.13. The nice instantaneous matching demonstrates the capability of the QS models in adapting to dierent wing kinematics, which is crucial for the application of the QS models on wing kinematics design. In Fig. 3.12(b), after 0.1 s, signicant discrepancies are observed in the horizontal force curves, which indicates the possible need of corrections for large wing root velocities. With the above results, the integrated simulation QS models show its capabilities in capturing the physics of the apping-wing MAV. Due to its computing eciency and iterative design capability, in the next subsection, the full dynamics model with QS models is quantitatively validated with the underwater tests. 77 Table 3.2: Parameter set for the QS model of the wing in underwater test. Name C L 1 C L 2 C L 3 C Dmax C D 0 C rot C rd C add C nc ^ x cp k h l h P 11 0.9 0.77 0.00 3.00 0.0 0.9 2.0 1.0 0.9 (2.21) 2.44 0.64 Table 3.3: Parameter set for the QS model of the body in underwater test. Name d e (mm) C b,d C b,a P 12 1.6 2.0 2.5 3.8.4 Experimental validation of the full dynamics model with underwater tests The results in the last subsections can partially re ect the accuracy and capability of the full dynamics simulation with QS models. However, the full dynamics model with CFD solver cannot provide more extra info on the inertial and conservative force terms than the QS models. To validate the entire full dynamics model presented in Sec. 3.33.6, the experimental tests were carried out under water with the robotic MAV II in Fig. 3.7(a). The experimental tests for the quantitative validation of the full dynamics model are presented in Figs. 3.14-3.15 and the 3rd-5th clips of the supplemental video S1.mp4. The CFD-calibrated QS parameters (with the integrated CFD simulation in Sec. 2.2.2) are listed in Table 3.2 and Table 3.3 for the wing and body, respec- tively. The parameters for the underwater modications are listed in Table 3.4. The geometric and inertial parameters for the wings and body are listed in Appendix B. For the free ight of the MAV in Fig. 3.14, the two wings ap symmetrically about thex b z b plane at 13 Hz with the nominal input prole described in Eqn. (3.57). Since there lacks a direct mapping from the electrical to the apping angle, the amplitude of the apping angle is estimated as 0 = 0:3 based on the initial wing behavior in the supplemental video S1.mp4. The bias parametersa 0 andb 0 are selected based on the experimentalz b prole. The gradient for the linear drifting isa 0 =1:0, and the 78 (b) Time (s) x b , z b (m) Planar coordinates of the flyer cenroid 0 0.5 1.0 1.5 2.0 2.5 0 0.04 0.08 0.12 0.16 x b -Sim x b -Exp z b -Exp z b -Sim θ b (rad) Pitch angle of the body 0.0 0.5 1.0 (c) Sim Exp (a) Sim Exp -0.5 Time (s) 0 0.5 1.0 1.5 2.0 2.5 Fig. 3.14. Full dynamics model validation with the underwater test for b 0 =0:05. (a) Dynamic behavior comparison between the experimental (left) and simulated (right) results. See the clip C4 in supplementary video S1.mp4 for more frames. (b) Instantaneous planar coordinates of the MAV body. (c) Instantaneous pitch angle of the body. 79 (b) x b , z b (m) Time (s) Planar coordinates of the flyer cenroid 0 1.0 1.5 2.0 2.5 0 0.04 0.08 0.12 0.16 0.5 x b -Sim x b -Exp z b -Exp z b -Sim (c) θ b (rad) Time (s) Pitch angle of the body 0 1.0 1.5 2.0 2.5 -0.5 0.0 0.5 1.0 Sim Exp 0.5 (a) Sim Exp Fig. 3.15. Full dynamics model validation with the underwater test for b 0 = 0:0. (a) Dynamic behavior comparison between the experimental and simulated results. See the clip C5 in supplementary video S1.mp4 for more frames. (b) Instantaneous planar coordinates of the MAV body. (c) Instantaneous pitch angle of the body. 80 Table 3.4: Parameter set for the model corrections in simulating underwater tests. Name v a limit (rad) k c (Nm) c c (Nm s) k leg (N=m) P 13 1.3 0.63 0.55 10 0.01 1.0 initial bias is b 0 =0:05. From Fig. 3.14, the simulated results highly resemble the experimental results. The matching is satisfactory, considering that only constant lumped QS and cor- rection parameters were applied to our simulation. To be specic, in Fig. 3.14(a), ve important video frames are plotted for each case. The full dynamics model is able to generate simulation results that capture the instantaneous physical behav- ior of the MAV. On the cycle-averaged behavior in Figs. 3.14(b)(c), the simulated results present the linear pattern of the vertical position, the parabolic pattern of the horizontal drifting, and the cycle-averaged prole of the body pitch angle. On the instantaneous behavior, the horizontal drifting is correctly simulated in terms of amplitude and phase. For the vertical position, the simulated results present stronger oscillations than the experimental prole. The reason possibly arises from the con- stant lumped correction, which cannot re ect local eects on velocity superposition. For the instantaneous pitch angle, the simulated prole shows a signicant discrep- ancy compared to the experimental test after t = 2.0 s. If the frames in Fig. 3.14(a) are checked, the experimental MAV has a non-zero yaw angle starting from t = 1.5 s. When the yaw angle is large, the perspective eect is strong enough to introduce errors to our 2D estimation of the body pitch angle. Though the parameter set for Fig. 3.14 provides us with reasonable predictions of the MAV dynamics, it is still possible that the results matching partially comes from the tuning of QS and correction parameters. To exemplify the calibrated set of parameters is capable of predicting new behavior of the MAV. Another experimental test was run with a nominal input with dierent e . Accordingly, the initial drifting bias is estimated as b 0 = 0:0. All the QS and correction parameters except b 0 are 81 the same as those for the rst test in Fig. 3.14. From Fig. 3.15, we can nd that the dynamics model can provide consistent accuracy for the body states when the apping bias is changed. For this test, the instantaneous results matching is even better than the case of b 0 =0:05. We have to point out that b 0 s are estimated based on the results in the video clips and the electrical signal sent to the actuators. The exact values are unknown without the identied signal mapping. Despite this uncertainty, the consistent accuracy on behavior prediction highly implies that our proposed full dynamics model is very eective in re ecting the dynamic physics of the MAV. At least, the experimental tests provide us with a preliminary validation for our full dynamics model. 3.8.5 Dynamic behavior analysis for MAV II ight in the air With the current design and fabrication quality of our robotic MAV II, the open-loop tests in the air cannot provide quantitatively repeatable results for model validation. The full dynamics model is still promising to be an eective tool for MAV's in-air behavior analysis. The results for three consecutive jump tests with = 110 Hz are shown in Fig. 3.16. These jump tests in the air were our validation tries be- fore the underwater tests. Here, we mainly emphasize the capability of our model in re-generating the behaviors of the MAV for dierent nominal inputs. The full dy- namics simulations were run with symmetric input signals for the two wings, and 2D trajectories are expected in the full dynamics simulation. For the baseline jump test, the nominal input signal has zero bias and is expected to y vertically to a high position. After an interactive process with the ground, the MAV jumps into a short height and then falls with a large body pitch angle from t = 0.46 s. The simulated results precisely re-produced the critical frames for this test. Then, to correct the body pitch angle, a large bias e was added into the second test. The new value of e over-corrected the pitch motion of the body into the other direction. The results captured the behavior of the MAV well before t = 0.07 s. However, after that frame, the MAV starts moving in a 3D trajectory with 82 0.0 s 0.22 s 0.10 s 0.46 s 0.72 s (a) 0.0 s 0.04 s 0.07 s 0.11 s 0.14 s (b) 0.0 s 0.06 s 0.12 s 0.19 s 0.25 s (c) Fig. 3.16. Comparisons of the experimental ying behaviors with the simulated results. (a) The baseline jump test with nominal bias b 0 = 0:0. (b) The test with large bias correction. (c) The test a small bias correction. 83 yaw rotations, which is not expected in our simulated results. In the third test, a small bias e was applied, and the simulated results are highly consistent with the experimental behavior inside the longitudinal plane. The above results re-emphasize our conclusion on the eectiveness of the full dynamics model on dynamic behavior prediction. In Fig. 3.16, the two-winged MAV presents a highly unstable pitch motion, which can be explained by the Eqn. (3.47). When there is a downward ow, R b provides a force component driving the body pitch further to destabilize itself, which is one reason for the inherent pitch instability of this conguration. 3.8.6 Eects of the added mass correction, initial conditions and hinge stiness To further demonstrate the capability of the full dynamics model and discuss uncer- tainties on the full dynamics model, three sets of simulations were implemented for the robotic MAV shown in Fig. 3.7(a). The simulations are for the free ight in the air with apping frequency 100 Hz and 0 = =3. To simplify the calibration process, these simulations used the QS parameters listed in Tables 3.1 and 3.3. The parameters for the pitching angle limit in Table 3.4 are also applied to simulations in this subsection. In Fig. 3.17(a), the results for two added mass correction c are shown. For the cases in the air, the in uence of the acceleration correction is very small for the overall dynamic response. The reason is that the acceleration from high-frequency apping is signicantly larger than that of the body motion. In Fig. 3.17(b), the ground interaction model (3.60a) is removed from the simulation. Compared to the case with ground model, the initial response of the no-ground-model one is signicantly dierent. However, similarities are still be observed from both the x b and b proles. In Fig. 3.17(c), multiple values of hinge stiness are tested in the simulation. In this simulation, the ground interaction model and hinge collision model are removed. 84 No ground model Ground model Vertical position of the body (a) (b) (c) Pitch angle of the body No ground model Ground model Pitch angle of the body -10 -8 Vertical position of the body α c = 0.63 α c = 1.00 Time (s) 0 0.05 0.10 0.15 0.2 0.25 0 x b (m) 0.02 0.04 -0.02 0 x b (m) 0.02 0.04 -0.02 Pitch angle of the body α c = 0.63 α c = 1.00 θ b (rad) -6 -4 -2 0 2 Vertical position of the body -0.10 -0.05 0 0.05 0.10 k h = 1.8 k h = 3.0 k h = 4.0 k h = 6.0 x b (m) Time (s) 0 0.05 0.10 0.15 0.2 0.25 Time (s) 0 0.05 0.10 0.15 0.2 0.25 k h = 1.8 k h = 3.0 k h = 4.0 k h = 6.0 Time (s) 0 0.05 0.10 0.15 0.2 0.25 Time (s) 0 0.05 0.10 0.15 0.2 0.25 Time (s) 0 0.05 0.10 0.15 0.2 0.25 θ b (rad) -6 -4 -2 0 2 θ b (rad) -6 -4 -2 0 Fig. 3.17. Eects of added mass correction, initial conditions and hinge stiness. (a) Sensitivity analysis on the added correction coecient. (b) Dynamic response with/without the ground interaction model for dierent initial conditions. (c) Dy- namic responses for dierent values of hinge stiness without the collision model and ground interaction model. For this robot, the simulation with k h achieved the highest vertical height. This exemplies an application of the formulated full dynamics model. 3.8.7 Towards new designs of apping-wing MAVs Behavior analysis and design applications of the full dynamics model for our two- winged conguration are shown in the last subsections. The study of the full dy- namics model also inspired us to investigate dierent congurations of apping-wing 85 (b) Thrust Lift Lift Drag (a) Fig. 3.18. New design of apping-wing MAVs inspired or potentially investigated by the studies of full dynamics models. (a) Photo sequence for the ight of Bee + , adapted from [18]. (b) Pitch oscillation of the body for low-frequency butter y-like ying. MAVs. In the procedure of exploring feasible experimental validation methods, the yaw damping and regulation were found not eective for both the open-loop and closed-loop tests. To conquer this issue, we rst improved the fabrication process to guarantee the symmetry of the MAV, which allows us to achieve the quasi-2D tests in a water tank in Sec. 3.8.4 and the high-quality open-loop jump tests in the air in Sec. 3.8.5. The ideas explored in the validation procedure were employed in the design of the yaw regulation for our four-winged robotic MAV, Bee + in Fig. 3.18(a) [18]. From the ight sequence, the quality of pitch motion and yaw motion were signicantly improved compared to the open-loop and closed-loop test clips for the two-winged conguration shown in the supplemental video S1.mp4. In our underwater tests, we employed the collision mechanism to limit the pitch angle of the wing. However, a hard collision between hinge shoulders for the MAV in the air was observed to induce strong initial yaw rotations, further unstable and complex ight trajectory. A new hinge design with soft collision to relieve this issue was proposed and will be presented in our other publications. For a butter- y MAV with a signicant lower apping frequency than 100 Hz (e.g., 20 Hz), the periodic 86 oscillation of the body is crucial for the lift generation for the MAV as illustrated in Fig. 3.18(b). With the instantaneous analyzing capability, the full dynamics model can be an eective tool for this kind of study. Therefore, no matter the physical insights from the validation process or the results from the design applications of our full dynamics model, these studies are to provide useful guidelines for the design of the new generation of the apping-wing MAVs. The high-quality experimental tests in Sec. 3.8.4 also have the potential to become benchmark data for other modeling research of apping-wing MAVs. 87 Chapter 4 Unied CFD Solver Prototype for Moving Boundary Problems In the previous chapters, dynamics models for a single wing and the entire MAV are simulated by integrating the equations of motion of the wing or MAV with the com- mercial Navier-Stokes solver, STAR-CCM+. Some drawbacks exist in the practical applications of this approach. First, the copyright license of a commercial package is not always available for all researchers in this eld. Second, as discussed in Sec. 3.8.1, massive parallel computing cannot be achieved for our integrated simulation with the commercial package. As the ODE solver for rigid body dynamics is very computationally ecient, the reason is attributed to the CFD solver. Also, a running of the simulation can frequently terminate because of the overset grid assembling fail- ure. Since we have no access to the source code of the package, the exact reasons for the low computing eciency or grid failure cannot be identied. Third, with a com- mercial package, the strongly coupled communication is not feasible. When data les are used for data exchange, the speed of communication and parallel implementation can slow down the entire computing work ow. Fourth, the rigid-body assumption for the wings is not always valid. From the setting-up procedure we know, it is not easy to implement a simulation in STAR-CCM+ with a exible wing actuated by a prescribed apping prole. These disadvantages prevent us from extending our 88 integrated simulation into a more ecient and generalized FSI framework. Great eorts were made by us to nd an alternative open-source solver to STAR- CCM+. Well-known open-source solvers with the overset grid capability include OVERFLOW [77], Overture CFD [78], NEK5000 [79] and the latest version of Open- FOAM [80]. OVERFLOW is a structured overset grid solver based on the high-order nite dierence method (FDM). As the patched grids present too many interfaces for complex geometry, OVERFLOW is not a good choice for our purpose. Over- ture CFD is the earliest CFD package with the overset grid capability. It solves the Navier-Stokes equations with the high-order FDM schemes for incompressible ow and second-order Godunov's scheme for compressible ow. During the use of Overture, its grid generator is not user-friendly, and the solver itself lacks interfaces for common mesh formats. Dierent from OVERFLOW and Overture, NEK5000 discretizes the Navier-Stokes equations with the spectral element method. The non-user-friendly in- terface also limits its application to common mesh formats. OpenFOAM is the most widely studied open-source CFD package. The overset grid approach is inherently integrated into its latest version, but still not robust. OpenFOAM denes an over- complicated and redundant input format, which is not friendly for quick setting-up. The supporting of standard mesh formats, such as CGNS, Fluent msh le, PLOT3D and VTK can signicantly improve the user-friendly property of the solver (e.g., SU2 [81]). The commonly used mesh generators (such as ICEM CFD, Pointwise and Gambit) for complex geometry provide output in these le formats. In terms of the capability for complex geometry and user-friendly interface, the commercial package such as STAR-CCM+ and Fluent outperform these open-source packages. The drawbacks of commercial and open-source packages and the need for a ex- ible FSI framework strongly motivated us to develop our in-house CFD solver for moving boundary problems. In this chapter, we rst discuss the strategies for the implementation of an overset grid-based CFD solver for moving boundary problems. Our research contributions on the overset grid assembling, automatic multigrid gen- 89 eration, and consistency enhancements for the segregated solver are addressed in Sec. 4.2 and 4.3. Testing examples for the CFD solver are nally presented to show its capability on enhancing algorithm accuracy and solving moving boundary problems. 4.1 Implementation strategies for the unied CFD solve prototype 4.1.1 Arbitrary Lagrangian-Eulerian approach and overset grid method FSI phenomena are inherently moving boundary problems (movable or deformable). For a control volume, the integral form of the Arbitrary Lagrangian-Eulerian (ALE) equations with a local moving frame of reference from [82, 83] is employed for the description of Navier-Stokes equations with moving boundaries, d dt Z d | {z } transient + Z @ (uu g )ndA | {z } convective = Z @ rndA | {z } diusive + Z @ q A ndA + Z q d | {z } source (4.1) where is the media density, u is the ow velocity, u g is the velocity of the surface element, is the transported scalar, is the dynamic viscosity, q A is the surface source vector, q is the volume source scalar, is the volume, A is the surface element and n is the local outward-pointing normal vector of the surface element. The symbol d dt denotes the non-local rate of change with location updating [82]. To be specic for our study, the incompressible, Newtonian and isothermal ow without gravity eect is considered and the discretized form for a specic mesh cell is d dt ( c;i ) | {z } transient + X A j j (u j u g;j )n j A j | {z } convective = X A j r j n j A j | {z } diusive X A j (p j i )n j A j | {z } diusive ; (4.2) 90 (a) (b) (c) (d) (e) Fig. 4.1. Common Lagrangian approaches for mesh boundary updating. (a) Im- mersed boundary method. (b) Dynamic mesh method. (c) Sliding grid method. (d) Overset grid method. (e) Boundary approximation with region domain mesh. where p j is the pressure of mesh face j and i is the Cartesian unit vector for . For the continuity equation, the right-hand side is zero. Schemes for the further discretization of Eqn. (4.2) are discussed in the next subsection. We rst discuss the updating strategies for the mesh boundary and associated unied implementation. In our research, the Lagrangian approaches for boundary tracking are preferred. The motion of the boundary is handled explicitly for local or global mesh updating [84]. Several common Lagrangian methods for updating moving boundaries are de- picted in Fig. 4.1. The immersed boundary method (IBM) was originally proposed for FSI simulation in biological ow problems [84, 85]. In the IBM method, the ow eld is typically discretized with the Cartesian grid, and the complex boundary is approximated by a non-conformal layer of mesh, as depicted in 4.1(a). Thus its coun- terpart for inviscid Euler equations is also referred to as the Cartesian grid method. For a moving boundary, a Lagrangian entity is treated as diuse [85] or as sharp 91 boundary with ghost cell or cut cell [84, 86]. For the dynamic mesh (or adaptive mesh) method [83, 87] in Fig. 4.1(b), the mesh is always conformal to the moving boundary and locally re-adjusted through morphing for a motion of the boundary in a small range. If the boundary moves in a large overall region, the global mesh must be re-generated, and the mesh can turn into highly distorted. Also, an interpolation algorithm is required for the transient term, and the accuracy can be further reduced. For a specic case with a spinning component, the sliding mesh in Fig. 4.1(c) can be employed with conformal interface geometry but non-conformal interface mesh. The sliding mesh is an advanced application of the interior interface concept available in most CFD packages [83, 88]. Fig. 4.1(d) is the overset grid method, which is the one for our integrated simulation. In the overset grid method (also called Chimera grid, composite grid, or overlapping grid method), the two or more mesh blocks are assembled together and present overlapping regions for data exchange. The overset grid method is very powerful for the representation of complex geometry and moving boundaries [82,89{92]. The overset grid method denes the interface relationship between two mesh blocks. The data structure can be directly used to express the conformal geomet- ric interface for the sliding mesh [93]. Since the dynamic mesh method is to update the grid of a single mesh block, it can be incorporated into the overset grid framework to re ect the local deformation of the boundary. The overall motion in a large range is expressed using the overset grid concept to avoid distorted mesh. In classic discus- sions of the immersed boundary method, the moving boundary is typically described as a continuous curve or surface [84{86]. If the mesh is imported from a general- purpose mesh generator, the boundary curve or surface is not easy to be obtained. If the boundary curve is replaced with a discretized mesh region as in Fig. 4.1(e), the hole-cutting process discussed in the next subsection for the overset grid method can be utilized in the approximation of the moving boundary for the IBM method. In our in-house solver development, the other three methods are incorporated 92 (a) (b) Hole-cutting Interpolation cell Host cell Donor cell (c) A B A Fig. 4.2. Overset grid assembly. (a) Original mesh blocks. (b) Articial boundary after hole-cutting. (c) Interpolation cell and stencil. into the implementation of the overset grid method. The application of the ALE equations and overset grid method (integrating the IBM, dynamic mesh and sliding mesh methods) form the \unied" concept for our in-house CFD solver with the capability of handling moving boundaries. 4.1.2 Overset grid assembling and discretization schemes for the segregated solver To discretize the ALE equations (4.1), the ow domain is rst discretized into a mesh block. If the overset grid technique is employed, the body-tting mesh block of the moving boundary and the background mesh block are dynamically assembled together as a mesh couple [82], as shown in Fig. 4.1(d) . The basic procedure for the overset grid assembling is described in Fig. 4.2. The background mesh A and moving mesh B blocks in Fig. 4.2(a) are prepared with a mesh generator. Part of the background mesh is removed from the region inside the outer interface of the moving mesh. An articial boundary is created as a layer of mesh colored in green in Fig. 4.2(b). The cells belonging to the articial boundary are referred to as the interpolation cells of the background mesh. For each interpolation cell in Fig. 4.2(c), its host cell from the moving mesh contains the centroid of the interpolation cell. The host cell and its neighbors constitute the interpolation stencil 93 (a) (b) c nb,j c i c i c nb,j c' i e j c' nb,j e' j e j n j Fig. 4.3. Finite volume discretization for a mesh cell. (a) Mesh neighboring infor- mation. (b) Gradient estimation with deferred correction, adapted from [95]. for the green dotted cell. For the moving mesh, its outermost layer of mesh cells is the interpolation-cell layer. As in our description, hole cutting and host search are the two essential steps of the overset grid assembling. The ecient host search algorithm from [82,94] and our proposed fast hole-cutting algorithm for unstructured grids are discussed in Sec. 4.2.1. The nite volume method (FVM) based on the collocated grid is employed for the discretization of the ow governing equations (4.2). The basic algorithms for the preliminary version of our in-house solver originate from [82, 83, 88, 95{97]. To solve discretized equations, the SIMPLE (Semi-Implicit Method for Pressure Linked Equations) [96, 98] algorithm is implemented to solve the momentum equations and continuity equation alternatively in a decoupled and iterative manner. Thus, our in-house FVM implementation is a segregated solver. In Fig. 4.3(a), c i is the current cell and its centroid, e j is the centroid of the face j, and c nb;j is the neighbor cell adjacent to face j and its centroid. n j is the outward-pointing normal vector . For the momentum equations, the transient term is discretized in the rst-order implicit Euler scheme, ( c;i ) n ( c;i ) n1 t =f n ; (4.3) where t is the time step size for the unsteady simulation and f is the abstract 94 function for spatial and source terms. The convective term in Eqn. 4.2 is discretized in the upwind scheme. For the face j of the cell c i j (u j u g;j ) n j A j = j (F j F g;j ); j = 8 < : c;i if (F j F g;j ) 0 nb;j if (F j F g;j )< 0 ; (4.4) whereF j is the velocity ux andF g;j is the grid motion ux. To satisfy the geometric conservation law, the grid motion ux is updated with the swept volume by the face [82,83,87,99] as F g;j = t ; (4.5) where is the volume swept by the face j in t time. If the mesh does not deform and the mesh block has no rotation, F g;j can be simply calculated as u g;j n j A j . The generalized procedure for the computation of is described in [82,99] For the diusive term, the second-order deferred correction method from [95] is utilized as, @ @n j = nb;j c;i (r nb;j r c;i ) n j + r old nb;j (r nb 0 ;j r nb;j )r old c;i (r c 0 ;i r c;i ) (r nb;j r c;i ) n j ; (4.6) where r is the position vector of the point in inertial frame andr old is the gradient of from last iteration. For the estimation of gradient at the cell centroid, the Gauss theorem is utilized, r c;i = P j j n j A j c;i : (4.7) We implemented the procedure described in the Sec. 3.3.1 of [82] to discretize the continuity equation or pressure equation. To complete the segregated solver, the contribution equations are derived for the interpolation cells with the stencil formed by the host and donor cells, as depicted in Fig. 4.2. The SIMPLE procedure 95 for the overset grid-based segregated solver with a strong coupling strategy for cell variables [82,88] is described as follows: 1. Discretize the momentum equations for the regular cells. a c;i c;i + N i X j=1 a nb;j nb;j =b c;i ; (4.8) where a c;i = t N i X j a nb;j and b c;i is the source term. 2. Interpolation equations of velocities for the interpolation cells. c;i + X k2I a k I;k = 0; (4.9) where I denes the set of the honor and donor cells. 3. Discretize the continuity equation as pressure correction. X j F 0 j = X j (F j F g;j );F 0 j = c;i a c;i j p 0 nb;j p 0 c;i d j n j jd j j; (4.10) wherep 0 is the pressure correction,F j is the face ux estimated from the corrected u j with u pred j . u pred j is the predicted velocity after solving the momentum equations (4.8) and (4.9). u j is dened in Eqn. (4.13). d j = r nb;j r c;i . 4. Interpolation cell updating for pressure correction with strong coupling strategy. p 0 c;i + X k2I a k p 0 I;k =p old c;i X k2I a k p old I;k ; (4.11) where p old is the pressure estimation from last iteration. 5. Correct pressure and velocity for the cells and faces p c;i =p old c;i + p p 0 c;i ; u c;i = u pred c;i c;i a c;i (rp 0 ) c;i : (4.12) where p is the under-relaxation correction. 96 Unstructured mesh .msh (Fluent) .cgns (NASA) Pointwise/Gridgen Gambit ICEM CFD In-house code Preprocessor .xyz (NASA) .vtk (Kitware) .bcs Sparse Matrix Lib + AMG solver Mesh Input files Boundary Metis + MPI Solver Output files .plt (Tecplot) .vtk (Kitware) Paraview Tecplot VTK library Postprocessor Overset grid toolkit Linux, C++ OgFSI Fig. 4.4. Technical strategies for the implementation of OgFSI. Our modied process for the segregated solver is discussed in 4.3. The discretiza- tion of our implementation is based on the collocated grids. To avoid the non-physical checker-board problem of the pressure eld [96], the Rhie-Chaw formulation is em- ployed by linking the face velocity with two pressure gradient estimations [82,95], u j = u pred j a c;i j p nb;j p c;i jdj j rp d j jd j j jd j jn j d j n j : (4.13) The incapability of this correction for small time step size and a revised version are discussed in Sec. 4.3. 4.1.3 Solver implementation and parallelization strategies The objective of our study in this section is to develop a general-purpose CFD solver for moving boundary problems (three-dimensional and incompressible laminar ow). The technical strategies are shown in Fig. 4.4. The solver is implemented in C++. The open-source library Metis is employed for partitioning the mesh. The formed linear system of equations from the overset grid FVM-based discretization is solved with an in-house sparse matrix library and a parallel algebraic multigrid (AMG) it- erative solver. The CFD solver is implemented based on the concept of unstructured mesh. The cell types include triangular, quadrilateral, tetrahedral, hexahedral, pyra- 97 Mesh files Partitionned mesh blocks Process 0 Process n Overset grid assembling FVM matrix Metis graph files Mesh files Partitionned mesh blocks Overset grid assembling Partitionned assembled blocks Parallel Multigrid generator FVM matrix t = t n Partitionned assembled blocks Parallel Multigrid generator AMG Solver AMG Solver Loop for i iter i iter > 0 Loop for t i iter > 0 i iter = 0 i iter = 0 Fig. 4.5. Parallel implementation with global indexing. mid, and mixed cell. The solver can read input mesh les in .msh, .cgns, .xyz and .vtk formats and output results in .vtk and .plt formats with in-house mesh import/export interfaces. The boundary condition le .bcs is a user-dened format. The provided interfaces to common mesh formats make the solver user-friendly for complex geom- etry and boundary conditions. The in-house overset grid toolkit was developed. The proposed fast overset grid assembling algorithm is discussed in Sec. 4.2. Our in-house solver is named as OgFSI (overset grid-based uid-structure interaction solver) with the capability for moving boundary and FSI problems. On the parallel implementation of the overset grid segregated solver, the global indexing strategy is employed as in Fig. 4.5. The strategy is for the multi-nodes and multi-process parallel computing with MPI. For each process, it reads in all mesh and partition les (Metis generated). Each process has a complete set of data storing all the cell information. The overset grid assembling operation is carried out for each process locally and independently. Thus, all the processes have the same global index for each cell and the overlapping cells. For each time stept n , when the automatically assembled grid blocks ready, the solver starts the parallel process. The parallel AMG generator only runs once for one time step at i iter = 0. Each process communicates with both its neighbor partitions on the same mesh block and the paired partitions from the other mesh for overlapping cells. The global indexing strategy can take large memory consumption, but each parti- tion can achieve fast access to the target cells. More importantly, with the global in- 98 (a) (b) n r d0 d P A B b out b in Q b h (c) Fig. 4.6. Distance-based hole cutting, neighbor-based linear host searching adapted from [82] and index re-numbering. dexing strategy, the diculty for overset grid assembling is signicantly reduced. The global indexing strategy permits a serial and easy overset grid assembling completed locally. This choice is a compromise between implementation diculty and comput- ing eciency. After achieving the capability of solving moving boundary problems, a parallel I/O le format and a parallel overset grid assembling algorithm are expected to be explored to solve the memory-eciency problem. 4.2 Fast overset grid assembling and automatic par- allel multigrid generation 4.2.1 Search-based fast overset grid assembling algorithm Distance-based hole cutting and neighbor-based host search For the overset grid assembling process, hole cutting and host search are the two crucial steps. When we were initializing our research, the distance-based hole cutting and host search algorithms from [82] were implemented as described in Fig. 4.6(a) and (b), respectively. In Fig. 4.6(a)he two mesh blocks are the background mesh A and the moving mesh B. Part of A is covered by block B. B has an internal boundary b in and outer boundary b out . b in is not used in the hole cutting procedure. The point 99 Table 4.1: Cell index re-numbering algorithm for each mesh block. 1. Import the cells with original indices. 2. Select a specic cell as the current cell, re-number the cell with the new index 0, set current level as 0, set the number of cells for current level as 1, and set n set = 0. 3. For the current level, traverse all the cells. (a). For the current cell, traverse all the neighbors of the current cell. If the neighbor cell is not assigned with a new index, assign the new index as n set and n set + +. The number of cells for next level also adds 1. (b). Set the next cell as the current cell. 4. Set the number of cells for next level as current level. Go back to 3 until n set equals the total number of cells. P represents the centroid of a cell that belongs to mesh A. To determine whether P is inside or outside b out , all faces on b out is traversed to nd the closest face centroid Q. The angle between face outward-pointing normal vector n and QP determines the location of P. The value d 0 denes the width the transition region with a ghost boundary b h . For A, all the cells inside b h are hidden from computation of the ow eld. The A cells closest to b h form the interpolation layer for A. The B cells closest to b out are the interpolation cells for B. When a cell is labeled as an interpolation cell (the red dot in Fig. 4.6(b), the host cell is searched with the trajectory depicted in Fig. 4.6(b) with the neighbor-based search algorithm from [82, 94]. Since the initial green cell can be randomly selected, if the connected line between the initial cell and the interpolation cell crosses b in or b h , part of the search trajectory may advance along the boundary cells. Cell index re-numbering If the mesh is generated with a general-purpose mesh generator, the indices for cells may not be well ordered as neighbors, as shown in Fig. 4.6(c). The automatic agglomeration procedure in the next subsection is a neighbor searching algorithm. If the neighbor cells have closer indices, the initial guess for the search algorithm can 100 be easier to be a good one, and the search procedure can be faster. A cell re-indexing algorithm is implemented as described in Table 4.1. To enable the re-numbering process feasible, the neighboring information (indices of the neighbor cells) must be available, which is the case for the Fluent mesh le. The re-numbering algorithm in Table 4.1 has a circular/spherical search trajectory from the center to the perimeter. For each cell, a limited number of searches are needed, and the searching algorithm is linear in terms of spatial complexity. Fast search-based hole cutting for large-scale problem In Sec. 4.2.1, the distance-based hole cutting and neighbor-based host search algo- rithms are described. Since the trajectory of the host search is always a polyline, it is a linear searching process for both 2D and 3D mesh blocks. It is computationally ecient for large-scale problems. However, for the hole-cutting process, if the mesh blocks are 2D with a small number of cells, the distance-based algorithm works pretty eciently, but when the mesh blocks are in high-resolution or a 3D simulation, the distance-based hole can in very low eciency. The reason is that, to determine a cell of A is inside or outside the outer boundary of B, all the face elements on b out must be traversed to nd the minimum-distance facet. The computing can be very slow for mesh blocks with million cells. A quick modication is to nd the minimum box containing the outer boundary b out to reduce the number of cells for computing, (a) (b) (c) (d) Fig. 4.7. Mesh topologies for overset grid assembling. (a) One layer with multiple moving objects. (b) Two layers. (c) Partially overlapping mesh blocks. (d) On mesh block partially outside background mesh. 101 d 0 n Q (a) (b) (c) (d) (e) (f) Fig. 4.8. Search-based fast hole-cutting algorithm.(a)-(c) Hole cutting for a bound- ary with continuous curvature. (d)-(f) Hole cutting for a boundary with non- continuous curvature. which is called the modied distanced-based hole cutting. For large-scale problems, this modication is still not fast enough. The neighbor-based host search and index re-indexing algorithms are linear search- ing algorithms. They are ecient enough for large-scale problems in our consideration. If the two algorithms can be combined, it is a very ecient hole-cutting algorithm. In our study, the overset grid approach is mainly employed for expressing moving boundaries instead of complex geometries, the simple topologies in Fig. 4.7(a), (b) are mainly considered. The more complicated mesh topologies in Fig. 4.7(c), (d). For the simple topologies, if the mesh resolution of the outer boundary is high enough compared to the mesh resolution of its background mesh, the fast search-based hole cutting can be achieved as in Fig. 4.8. Algorithm steps are described in Table 4.2. To test the eciency of the fast hole-cutting algorithm, we run the assembly algorithm for the 3D mesh in Sec. 4.4 with around 1.5 million cells for the background mesh and 1.4 million cells for the wing mesh, the speed testing for one-step assembling 102 Table 4.2: Proposed fast search-based hole-cutting algorithm. 1. For each face element on b out , its centroid Q can be push inward with the normal vector n and transition distanced 0 to the green dotted position as in Fig. 4.8(a). 2. The neighbor-based linear search is utilized to nd the host cell of the green dotted point from mesh block A as in Fig. 4.8(a). 3. The resolution of b out is high enough, at least one complete layer of mesh in red can be found out to isolate the outer and inner of the dashed articial boundary as in Fig. 4.8(c). If a seed cell adjacent to the red layer is inside the inner of the articial boundary, for most cases, the index re-numbering algorithm can be applied to mark all the cells as hidden cells. 4. For a special case with non-continuous gradients for a boundary as in Fig. 4.8(d), the blue cell cannot indicate the inner or outer region of the articial boundary properly. Wrong hole-cutting is possible to happen. If the mesh block is solid mesh without inner boundary, the centroid of the mesh block can be used to nd the seed cell. If the block is moving boundary with an inner wall, the seed cell can be found as a cell adjacent to the boundary as in Fig. 4.8(e). 5. After the seed cell is identied correctly, run the index re-numbering algorithm to mark all hidden cells. 6. Traverse all the red cells. If the cell is not adjacent to a hidden cell, it is marked as a regular cell. The identied articial boundary is as in Fig. 4.8(f). on the same platform (Linux and Intel i7-9700K) is listed in Table 4.3. The search- based algorithm is signicantly faster than the distance-based algorithm. Though not feasible to time the non-initializing assembly process in STAR-CCM+, our proposed search-based algorithm is at least competitive to the commercial code. For the applications of the fast search-based hole-cutting algorithm to other cases, if the mesh resolution of the outer boundary is not high enough, the dashed articial boundary can be found through a line or plane interpolation for auxiliary points in 103 Table 4.3: Eciency testing for the hole-cutting algorithms. Distance-based Modied distance-based Search-based STAR-CCM+ 80 s 10 s 0.1 s 10 s (initialize) < 1 s (re-assembling) space. For the moving boundary problems, the body-tting mesh typically owns a higher resolution than the background mesh. 4.2.2 Hybrid parallel automatic multigrid generation For the linear system formed from the discretization of Navier-Stokes equations using the nite volume method, the equations are assumed in the form of a c;i c;i + N i X j=1 a nb;j nb;j =b c;i : (4.14) Here, Eqn. (4.8) is re-presented for quick access. For the row from the descritization of overset grids, the contributions equations are in the same form with the diagonal coecient as \1", i.e., Eqn. (4.9). There are some associated special topology features and numerical diculties including: A n-1 x n-1 = b n-1 A n-2 x n-2 = b n-2 Algebraic agglomeration (a) (b) Fig. 4.9. Matrix topology and AMG process for the overset grid-based solver. 104 P0. From Fig. 4.9(a), the matrix from ow eld discretization is not symmetric due to the interpolation cell rows (the green boxes). Some special fast algorithms such as conjugate gradient method cannot be used directly. P1. The hidden cells are nominally not involved in the ow eld computing [82]. However, to simplify the matrix formulation, the rows for the hidden cells (the red boxes in Fig. 4.9(a) still exist in the computing matrices with the only non-zero entry 1. No neighbor cell information is available for these rows. P2. For the interpolation cells (the green boxes in Fig. 4.9(a)), the non-zero non- diagonal entries are interpolation coecients across two mesh blocks. No neighbor cell's information can be extracted from the non-zero entries. P3. If two moving mesh blocks partially overlap each other as the ones Fig. 4.7(c), the matrix topology can change signicantly due to the movement of meshes. The formulation of the multigrid can become very complicated. To simplify the discussion of the automatic multigrid generation, we limit the use of the complicated case relating to P3. In an assembly of overset grids, one upper layer mesh block can always be fully contained by the paired background mesh block. Automatic agglomeration algorithm and serial generator In the algebraic multigrid procedure, a matrix for a coarser layer is formed automat- ically from the matrix of the upper ner layer as depicted in 4.9(b). To simplify the notation, for each layer of the multigrid, a nb;j represents the non-zeros entries at the non-diagonal positions. For the entries corresponding to the same mesh block, the most straightforward agglomeration condition for determining the neighbor relations is selected as a nb;j >: (4.15) If the condition (4.15) is satised, the cellc nb;j is treated as the neighbor of c i . For the matrix rows relating to the diculties P1 and P2, there is no neighbor information 105 available in the matrix itself. Extra geometric info must be passed from the geometric mesh to the multigrid generator as shown in Fig. 4.10(a). Otherwise, for P1, each hidden cell becomes an independent cell in all coarser meshes and for P2, improper neighbor relations are possibly applied in the coarser level of meshes. To solve P1, using the extra geometric info, two methods can be considered to implement: (1). Remove the hidden regions from the computing matrix to construct a set of reduced matrices. This can reduce the dimensions of the computing matrix. However, for moving boundary problems, the indices between the cell index and the matrix rows number must be updated at each time step. Due to this mapping, the parallelization requires more information to be shared and communication becomes less ecient. (2). Since the overlapping information can be transferred to the multigrid constructor, all the hidden cells can be united as one single cell for each coarser layers as in Fig. 4.10(b). Less storing memory is required. With this method, the cell index and the corresponding row index is the same. The communication interfaces between partitions are similar to that of the non-overset grids. For P2, the hybrid geometric-algebraic procedure is employed by passing the list indicating the overlapping status to the multigrid solver. In this procedure, the inter- A s AMG generator Unifying (a) (b) (c) (d) Fig. 4.10. Automatic agglomeration algorithm for serial and parallel generators. 106 Table 4.4: Serial search algorithm for the coarser layer of grid. 0. Store the overlapping status of each cell in a 1D array, and pass this array to the multigrid generator for the coarser layer with the discretized matrix. 1. Traverse all the matrix rows (cells) in each mesh block and the default group size of a coarser cell is 4. (a). If current cell is a regular cell and one neighbor cell is already added into a non-full coarser cell. Add the current cell into this coarser cell. (b). If current cell is a regular cell and no neighbor cell is added into a coarser cell. Create a new coarser cell. Add its neighbor cells into the new coarser cell until it is full. During this step, if the neighbor cell is an interpolation cell and the new coarser cell is full, adjust the group size larger than 4. (c). If an interpolation cell is added into current coarser cell. Mark the current coarser cell as an interpolation cell and store the overlapping status. 2. Traverse all the matrix rows (cells) in each mesh block again. (a). If the neighbor cell is an un agged interpolation cell and in the same mesh block with current cell, add the interpolation cell into corresponding coarser cell. 3.Traverse all the matrix rows (cells) in each mesh block again. (a). If one cell is not agged, create an independent coarser cell. (b). Unify all hidden cells as one coarser cell for each mesh block. polation cell cannot be the rst cell for agglomeration, and it can only be absorbed by an existing coarser cell. Also, each cell (matrix row) can only be agglomerated into the same mesh block. No inter-block coarsening is allowed. Thus the automatic agglomeration procedure in Table 4.4 is feasible. The auxiliary information is in Fig. 4.10(c). In Table 4.4, the matrix row is equivalent to the concept of a cell. The procedure is to form the computing matrices for the coarser layer. Parallel agglomeration searching In the parallel version of the solver based on domain decomposition, for each partition, there are two types of neighboring partitions: the adjacent neighbor partition from the 107 Table 4.5: Parallel search algorithm for the coarser layer of grid. 0. Store the overlapping status of each cell in a 1D array and pass this list to the multigrid generator for the coarser layers generation with the discretized matrix. 1. Traverse all the matrix rows (cells) in each partition and the default group size of a coarser cell is 4. (a). If the current cell is a regular cell and one neighbor cell is already added into a non-full coarser cell. Add the current cell into this coarser cell. (b). If the current cell is a regular cell and no neighbor cell is added into a coarser cell. Create a new coarser cell. Add its neighbor cells into the new coarser cell until it is full. During this step, if the neighbor cell is an interpolation cell and the new coarser cell is full, adjust the group size larger than 4. (c). If the interpolation cell is added into current coarser cell. Mark the current coarser cell as an interpolation cell and store the overlapping status. 2. Traverse all the matrix rows (cells) in each partition again. (a). If the neighbor cell on the same partition is an interpolation cell and is not agged, add the interpolation cell into the corresponding coarser cell. 3. Traverse all the matrix rows (cells) in each partition again. (a). If one cell is not agged, create an independent coarser cell. (b). Unify all hidden cell as one coarser cell for each partition. 4. Communicate between partitions and mesh blocks to construct the global indexing. same block and the interpolation neighbor partition from the other mesh block. To implement the automatic agglomeration procedure in Table 4.4, the cells are coarsened locally on each partition and inside their own mesh blocks. The coarsening data information is constructed for both sets of neighbors separately. The parallel search algorithm is listed in Table 4.5. The auxiliary information is in Fig. 4.10(d). In comparison to the serial version, one more step is added to handle the communication to construct global indexing. Specially, step 3 is to guarantee the correct coarsening of the interpolation cell which is isolated by hidden region and partition boundaries. 108 4.3 Smoothness and consistency enhancements OgFSI was rst implemented with the algorithms presented in Sec. 4.1. This imple- mentation is an old version for results comparison. Thus it is referred to as OgFSI-old later. Numerical tests were implemented with the lid-driven cavity, channel ow and cylinder ow problems. Specic issues related to the regularity of the mesh, the time step size and the overset grid coupling strategy are discussed with these examples. Revised schemes are proposed to resolve or interpret the mentioned numerical issues. Corresponding improved results are presented in Sec. 4.4. 4.3.1 High-order face interpolation for non-regular mesh To test the numerical consistency of OgFSI-old on non-regular grids, the unsteady channel ow problem with partially distorted mesh was run as in Fig. 4.11(a). The time step size is t = 0:1 s. The boundaries include a velocity inlet (u in = 0:01 m/s, h = 0:3 m, l = 2 m), a pressure outlet and no-slip upper and bottom walls. Non- physical discontinuities are observed in Fig. 4.11(b) in the middle of the channel. To understand these results, the same settings were also run with STAR-CCM+ and Fluent. Similar non-smooth features of the velocity eld are observed in Fig. 4.11 (c) and (d). The reason was attributed to the non-regularity of the mesh. In [100], the non-regularity issue of the mesh is discussed. The non-smoothness was identied to be the updating scheme of face variables after cell variables are corrected in Eqn. (4.12). In the implementation of OgFSI-old, the rst-order scheme was employed as u e = 0:5(u 1 + u 2 );p e = 0:5(p 1 +p 2 ); (4.16) where the denitions of variables are in Fig. 4.12. The updated variables with Eqn. (4.16) are not the values for the face centroid e. The high-order scheme with deferred 109 (a) (b) (c) (d) (e) h l Fig. 4.11. Non-smoothness of velocity eld for non-regular mesh. (a) Non-regular mesh for channel ow. (b) Result of velocity eld from OgFSI-old. (c) Result from STAR-CCM+ 12.02. (d) Result from Fluent 15.2. (e) Result from OgFSI with a high-order interpolation scheme. correction is utilized to compensate the location error, u e;k = (1)u 1;k +u 2;k +ru old e 0 ;k r e 0 e ; (4.17a) p e = (1)p 1 +p 2 +rp old e 0 r e 0 e ; (4.17b) where k denes the coordinate direction, and = r 1 r d r d r d . Improved result for the channel ow using this high-order scheme is presented in Fig. 4.11(e). 110 u 1 , p 1 u 2 , p 2 e e' r d r 1 r 2 r e'e u e , p e Fig. 4.12. High-order interpolation for the updating of face variables considering the non-regular shape of the cell. 4.3.2 Smoothness and consistency of FVM for small t In Sec. 2.2.4 and 3.8.1, integrated simulations for a wing apping at 100 Hz were run with a time step size 0.1 ms. OgFSI-old was tested with a xed-cylinder ow (d = 0:1 m). Structured multi-block mesh in Fig. 4.13(a) is employed for the discretization of the ow eld. Results of the pressure eld at the 2nd time step are presented for two dierent ts. For t = 0:1 s, smooth pressure eld is achieved. In the contrast, highly non-smooth pressure distribution exists for t = 0:1 ms. The same simulation setup at t = 0:1 ms was run with STAR-CCM+ and Fluent. Worse result was observed in STAR-CCM+ because when more time steps were run and the simulation diverged. Fluent was able to achieve smooth pressure eld for the small t = 0:1 ms. In Fig. 4.13(c), non-smooth feature in region A resembles the non-smoothness in Fig. 4.11(b). High-order schemes were tried for the discretization of velocity, pressure, face values, but failed. When the non-smoothness in region B was checked, the feature is highly similar to that of the \checker-board" problem. In [95], it is pointed out that the pressure correction can be ineective for small t with the Rhie-Chow interpolation. The Rhie-Chow formula is restated as below, u j = u pred j a c;i j p nb;j p c;i jdj j rp j d j jd j j ! jd j jn j d j n j ; (4.18a) a c;i = t N i X j a nb;j : (4.18b) 111 (a) (b) (c) (d) (e) A B (f) d Fig. 4.13. Non-smoothness of the pressure eld at the 2nd time step. (a) Multi-block mesh for the cylinder ow. (b) OgFSI-old at t=0.1 s. (c) OgFSI-old at t=0.1 ms. (d) STAR-CCM+ at t=0.1 ms. (e) Fluent at t=0.1 ms. (f) OgFSI at t=0.1 ms. If t is in a very small value,a c;i is in a large value. The 2nd term of the right-handed side of (4.18a) for pressure gradient correction can be ineective. The suggested scaling-up correction with a constant number from [95] was tried, but did not help. The dependence of the Rhie-Chow correction on the time step size was discussed in [101{104], but the improved Rhie-Chow interpolation schemes were not able to solve the presented non-smoothness problem in Fig. 4.13(c). 112 The \checker-board" problem is a spatial non-smoothness issue. Since the pressure eld in the 1st time step is smooth, we tried to remove the transient terms from the predicted velocity u pred j as u pred,s j = u pred j u n1 j ; (4.19) where u n1 j is the face value from last time step. The new predicted face velocity is the revised one for the updating of F j in Eqn. (4.10). Using this correction, smooth results are achieved as plotted in Sec. 4.4. Though not rigorously proved, the non-smoothness issue can be explained with quantitative analysis on the error propagation of the continuity equation (cell centroid velocity) or conservation of mass (face centroid velocity). For the time step n, the continuity equation error residual from scheme (4.19) is @u n c;i @x + @v n c;i @y =R 1 + n X k=2 R k : (4.20) For the smooth simulation results presented in 4.13(f), the step residualR 1 is around 0.1 and R k for k > 2 is around 1e-4 with its specic units. Also, the sum of R k is not observed to increase in a long-time running. The errors from the 1st time step dominate the error residuals for the continuity equation or pressure correction equation. We run the same simulation in Fluent for a case with long physical time, and the error residual of the continuity is also close to being a constant value (around 0.9), which is equal to that of the 1st time step. The results are consistent with those from OgFSI, which indicates that the pressure eld is not able to be corrected properly if the velocity errors from the 1st time step are propagated to later time steps. To achieve consistency of pressure correction for small time steps, the pressure eld in the 1st time step must be solved accurately with small absolute error residuals. 113 4.3.3 Coupling strategies for interpolation cell of the overset grids (a) (b) (c) (d) Fig. 4.14. Non-smoothness of the velocity and pressure elds for overset grid solu- tions at small time step t = 0:1 ms. (a) Velocity eld by Fluent. (b) Pressure eld by STAR-CCM+. (c) Smooth velocity eld by OgFSI. (d) Smooth pressure eld by OgFSI. After using the segregated overset-grid process described in Sec. 4.1.2, non-smooth features can be observed in the velocity and pressure elds for the cylinder ow. The mesh topology for the overlapping region is similar to that in Fig. 2.9. To illustrate the generalities of this issue, the velocity eld from Fluent (coupled solver) and the pressure eld from STAR-CCM+ (segregated solver) are plotted in Fig. 4.14(a) and (b), respectively. Non-physical spiking regions appear in both gures, which is consistent with the results from OgFSI-old. To interpret these non-smoothness problems, dierent interpolating schemes for the overlapping regions were tried. The 114 (a) Interpolation face Host cell Donor cell (b) Interpolation cell Host cell Donor cell Fig. 4.15. Strong coupling vs loose coupling for solving pressure correction. (a) is re-plotted from Fig. 4.2(c). non-smooth features were identied to arise from the strong coupling of the data exchange among mesh blocks. To resolve the above non-smoothness problems, as shown in Fig. 4.15(b), a loose coupling pressure correction scheme for the interpolation cells is employed as a p 0p 0 c;i + N J X k2J a k p 0 J k + N I X k2I a k N W X l2W l p 0 W l =b p 0; (4.21) where J is the set of faces which are adjacent to regular cells, I is the set of the interpolation faces, andW is the set of donors cells for faces inI. The loose coupling interpolation scheme, Eqn. (4.21), is to replace Eqn. (4.11) in the step 4 of the segregated SIMPLE procedure in Sec. 4.1.2. As depicted in Fig. 4.15(b), the two mesh blocks exchange the information of pressure correction through face values instead of cell values, which is the source for the \loose coupling" concept. Smooth velocity and pressure ow elds using OgFSI with the loose-coupling interpolation scheme (4.21) are plotted in Fig. 4.14(c)(d). 4.4 Testing examples for OgFSI In this section, numerical examples are presented for OgFSI to show its capability on assembling overset grids, automatically constructing multigrid layers, and simulating unsteady and moving boundary problems. 115 (a) (b) (c) (d) (e) (f) Fig. 4.16. Overset grid assembly with the search-based fast algorithm. 4.4.1 Overset grid assemblies with search-based algorithm The overset grid assembling algorithm in Sec. 4.2.1 was implemented for several dierent mesh types. Fig. 4.16(a) is the one-layer topology with one xed back- ground mesh and one moving cylindrical mesh. Fig. 4.16(b) is the example of a solid rectangular moving mesh assembled with one background mesh. Fig. 4.16(c) is the two-layer topology that resembles the mesh topology for the full dynamics simulation in Sec. 3.5. Fig. 4.16(d) is the one-layer mesh topology with two moving bodies. Fig. 4.16(e) is an example of parallel assembled grids. Fig. 4.16(f) is an example of 3D assembled grids for a wing plate. 4.4.2 Cell index re-numbering and multigrid generation Fig. 4.17 shows two examples of cell index re-numbering and automatic multigrid coarsening for single-block grids. The color represents the index of each cell. Blue 116 (a) (b) Fig. 4.17. Cell index re-numbering and coarsening for single block mesh. color represents a small value and red color represents a large value. Cells in the same color belong to the same coarser cell. In Fig. 4.17(a), an unstructured grid is gen- erated using ANSYS ICEM CFD with original indices almost randomly distributed. To re-number cell indices, the new zero cell is selected around the center of the mesh block, and the circular coarsening pattern is observed. Fig. 4.17(b) is generated for the structured grid and the original indices are ordered vertically. The cell with the new zero index is at the bottom left corner. The cells are coarsened diagonally. Fig. 4.18 presents the results of the automatic agglomeration process for the assembled overset grids. The cell index for each mesh block starts from zero indepen- dently. The rst image in Fig. 4.18(a) shows the original mesh structures with blue lines. In the assembled grids, the overlapping region between the two mesh blocks can be observed. For each block, the mesh is coarsened independently. The rst image in Fig. 4.18(b) shows the decomposed partitions with the METIS library. In the parallel case, all the mesh cells are coarsened locally on each partition. Fig. 4.18(c) and (d) demonstrate the capability of our multigrid generator for the cases with multiple mesh blocks and 3D mesh blocks. The last image in Fig. 4.18(c) shows the geometric conforming of the coarser grids. The last two images of Fig. 4.18(d) shows the internal structure of the 3D coarser grids. From the above examples, the proposed automatic multigrid generation algorithm 117 (a) (b) (c) (d) Fig. 4.18. Examples of the automatic multigrid generator for overset grids. (a)- (b) Automatic agglomeration for serial and parallel coarsening. (c)-(d) Multigrid generation for multiple mesh blocks and 3D meshes. is demonstrated to be very eective for the unstructured parallel overset-grid solver, which is one basis for the ecient parallel AMG solver for OgFSI. 4.4.3 Finite volume solver based on overset grid The accuracy of OgFSI is rst validated for the single-block mesh with the lid-driven cavity problem. Interpolation steps in Sec. 4.1.2 are neglected in the implementation. Dimensions for the geometry are 1.0 m 1.0 m and the driven velocity is 0.001 m/s. The operating Reynolds number isRe 68:5. The steady results for mesh resolution 118 OgFSI STAR-CCM+ Fluent A A A OgFSI, triangular A A A A (a) (c) OgFSI STAR-CCM+ OgFSI (IBM) (b) (d) Fig. 4.19. x velocity eld for the lid-driven cavity problem. (a) Single structured quadrilateral mesh block with 100 100 cells. (b) Single unstructured triangular mesh block. (c) x velocity eld and mesh assembling for OgFSI and STAR-CCM+ with overset grids. (d) FVM-based IBM method using the overset grid hole-cutting. 100 100 from OgFSI, STAR-CCM+ and Fluent are plotted in Fig. 4.19(a). From the visualized x velocity elds, all solvers generate very close results and the minor dierences exist in region A. To illustrate the capability of OgFSI on solving ow problems with dierent mesh types, the results with unstructured triangular mesh are plotted in Fig. 4.20(b). Since the results are not smooth enough, the x velocity eld shows the traces of the triangular mesh. Results on the steady problem demonstrate the correctness and accuracy of the spatial discretization. The results for the stationary overset grid and immersed boundary methods are plotted in Fig. 4.19(c) and (d), respectively. An additional xed wall is added at the 119 center of the lid-driven cavity. The diameter of the wall is 0:2 m and the external diameter of the body-tting mesh is 0:5 m. The x velocity and mesh assembling agreements between OgFSI and STAR-CCM+ illustrate the proper implementation of OgFSI. However, the two solvers show some dierences in region A where OgFSI generates a slightly larger low-velocity area. Though both solvers employ the rst- order distance weighted formulas for the interpolation cells, the specic interpolation t = 0 s t = 50 s t = 100 s t = 150 s (a) (b) t = 0 s t = 10 s t = 20 s t = 30 s t = 40 s (c) (d) Fig. 4.20. Results for the moving rigid overset grid solvers for unsteady simulation in 150 s. (a) Rigid cylinder with STAR-CCM+. (b) Rigid cylinder with OgFSI. (c) Deformable cylinder with OgFSI. (d) Mesh deformation for the cylinder with OgFSI. 120 stencil for each cell can be dierent. The width of the transition region is also dierent. We attribute the discrepancies between the two solvers to the actual realization of the overset grid algorithm. The results of the immersed boundary based on the overset grid hole cutting provides evidence for the generality of our unied framework idea. At this Reynolds number, the non-conformal feature of the immersed boundary does not prevent the generation of a reasonable x velocity eld. The validation of OgFSI for unsteady problems is presented together with the moving boundary results in the next subsection. 4.4.4 Results for moving cylinder problems To validate OgFSI for moving boundary problems, the results for rigid cylinders are rst presented in Fig. 4.20(a)(b). The cylindrical wall is moving in a prole u = 0:0025 cos(0:0126t) m/s. OgFSI generates a smooth velocity eld that matches well that from STAR-CCM+. More interestingly, if we compare the rst frame (t = 0 s) between Figs. 4.20(a) and (b), OgFSI has non-physical spikes that appear in the results of STAR-CCM+. Figs. 4.20(c)(d) present results for deformable moving boundary problems. The moving cylindrical wall is allowed for local deforming ind x (t) =d(t 0 )0:04 sin(2ft). x velocity eld and the corresponding deformable mesh are plotted in Fig. 4.20(c) and (d). Here, the deformation of the mesh is prescribed as a function of time. If a structure solver is included, the deformation shall be solved dynamically with the equations of elasticity. 4.4.5 Results for moving wings To illustrate the capability of OgFSI on simulating apping-wing problems, the mov- ing wings are plotted in Fig. 4.21. OgFSI is able to simulate both the rigid and deformable wing problems. 121 (a) (b) (c) (d) (e) (f) Fig. 4.21. x velocity elds for the moving wings with OgFSI. (a) Rigid 2D elliptic wing with prescribed horizontal velocity and pitching angular velocity. (b) Rigid 3D rectangular wing. (c)-(d) Vertical elliptic wing in horizontal motion with enlarged minor axis. (e)-(f) Moving elliptic wing with local bending. 122 Chapter 5 Conclusions and Future Work Conclusions In this dissertation, a set of multi-delity dynamic models and tools were proposed and implemented for the apping-wing MAVs. The results proved that these models are able to be useful tools for the aerodynamic, mechanical and controller design of the apping-wing MAVs. The development of OgFSI provided an overset grid- based solver prototype with consistency improved algorithms for the moving boundary problems. We can achieve more exibilities in the dynamic simulations for apping wings and other FSI investigations. The integrated FSI simulation for a single wing was rst proposed to estimate aerodynamic forces on a single wing in the passive pitching motion. In the integrated simulation, the equation of motion and the Navier-Stokes equations are solved sep- arately and sequentially following a weakly coupled FSI assumption. To achieve a fast estimation of the wing dynamics, quasi-steady models were incorporated into the equations of motion. Revised formulas were proposed to improve the adaptability of the QS models to dierent wing kinematics. The pitching angle limit placed by hinge geometry is preliminarily modeled using a geometric limiter or an additional spring- damper subsystem. The underwater simulation for a wing apping at low frequency 123 exemplied the condition for the weak coupling assumption. An alternative form of the Lagrange's equations was derived by removing repeated terms from the inertial force expressions. The new compact form allows us to man- ually derive equations of motion for complex system with high numerical eciencies. The full dynamics model of the entire MAV was rst integrated with a proprietary CFD solver in a quasi-high order ODE process. The demand on high computing e- ciency motivated us to extend the quasi-steady models for a single wing to the entire MAV by considering the wing root kinematics. The quasi-steady models were rst numerically validated with CFD results. The full dynamics model with quasi-steady formulas was then experimentally validated with quantitative underwater tests and qualitative jump tests in the air. The validation results showed the eectiveness of the full dynamics model on predicting the dynamic behavior of apping-wing robotic MAVs. To achieve the maximum exibility and accessibility of CFD computing, the in- house unied CFD solver (OgFSI) based on the overset grid method was developed for FSI and moving boundary problems. The fast search-based hole-cutting algorithm, the hybrid parallel multigrid generator, and the smoothness improved discretization enabled the segregated solver to be ecient and eective in tackling moving boundary problems. Future work In modeling a single wing and the entire MAV, the mapping between the electrical input signal and the apping angle of the wing is assumed to be in linear time- invariant model. The apping prole is the same as that of the electrical signal prole. However, the underwater experimental results showed a strong interaction between the actuator-transmission and the water-wing system. The apping prole is challenging to be linearly mapped from the input electrical signal. To simulate 124 the wing dynamics more independently without observation from experiments, the actuator-transmission model must be part of the wing dynamics. From the design perspective, the actuation-transmission model is also useful for the mechanical design of the actuator for the robotic MAV in the air. In this dissertation, most results on the full dynamics model were presented for the purpose of quantitative or qualitative validations. Intuitively, the next step work would be to explore the applications of the full dynamics model in mechanical and con- troller design of the apping-wing MAVs. Also, the apping-wing MAV is a nonlinear multibody system, whose full dynamics model can serve as a platform for nonlinear dynamic phenomena investigation. The applications of OgFSI on the dynamics of a single wing and entire MAV are expected in my post-Ph.D. research. The development of a general purpose CFD solver is the ultimate long-term objective. Some promising functionalities to be added into OgFSI include memory-ecient algorithms for parallel overset grid assembling, turbulence models for moderate and high Reynolds numbers, and the integration of nite dierence and nite element methods into this solver. 125 Bibliography [1] T. Weis-Fogh and M. Jensen. Biology and Physics of Locust Flight. I. Basic Principles in Insect Flight. A Critical Review. Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences, 239(667):415{458, 1956. [2] T. Weis-Fogh. Quick Estimates of Flight Fitness in Hovering Animals, Includ- ing Novel Mechanisms for Lift Production. Journal of Experimental Biology, 59(1):169{230, 1973. [3] C. P. Ellington, C. Van Den Berg, A. P. Willmott, and A. L.R. Thomas. Leading-Edge Vortices in Insect Flight. Nature, 384(6610):626, 1996. [4] S. P. Sane. The aerodynamics of insect ight. Journal of Experimental Biology, 206(23):4191{4208, 2003. [5] T. Weis-Fogh. Biology and Physics of Locust Flight II. Flight Performance Of the Desert Locust (Schistocerca gregaria). Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences, 239(667):459{510, 1956. [6] M. H. Dickinson, F. O. Lehmann, and S. P. Sane. Wing Rotation and the Aerodynamic Basis of Insect Flight. Science, 284(5422):1954{1960, 1999. [7] J. P. Whitney and R. J. Wood. Aeromechanics of Passive Rotation in Flapping Flight. Journal of Fluid Mechanics, 660:197{220, 2010. [8] E. K. Singer, L. Chang, A. Calderon, and N. O. Perez-Arancibia. Clip-Brazing for The Design and Fabrication of Micro-Newton-Resolution, Millimeter-Scale Force Sensors. Smart Materials and Structures, 2018. [9] G. K. Taylor and A. L. R. Thomas. Dynamic Flight Stability in the Desert Locust Schistocerca gregaria. Journal of Experimental Biology, 206(16):2803{ 2829, 2003. [10] M. Sun and Y. Xiong. Dynamic Flight Stability of a Hovering Bumblebee. Journal of Experimental Biology, 208(3):447{459, 2005. 126 [11] C. T. Orlowski and A. R. Girard. Modeling and Simulation of Nonlinear Dy- namics of Flapping Wing Micro Air Vehicles. AIAA Journal, 49(5):969{981, 2011. [12] S. Chang and Z. J. Wang. Predicting Fruit Fly's Sensing Rate with Insect Flight Simulations. Proceedings of the National Academy of Sciences, 111(31):11246{ 11251, 2014. [13] S. Deng, M. Percin, B. W. van Oudheusden, A. Tenaglia, and B.D.W. Remes. Force and Flow Field Measurements of a Bio-Inspired Flapping-Wing MAV `DelFly Micro' in Hovering Flight. In 32nd AIAA Applied Aerodynamics Con- ference, pages 16{20, 2014. [14] M. Kar asek, F. T. Muijres, C. De Wagter, B. D.W. Remes, and G. C.H.E. de Croon. A Tailless Aerial Robotic Flapper Reveals that Flies Use Torque Coupling in Rapid Banked Turns. Science, 361(6407):1089{1094, 2018. [15] A. Roshanbin and A. Preumont. Yaw Control Torque Generation for A Hovering Robotic Hummingbird. International Journal of Advanced Robotic Systems, 16(1):1729881418823968, 2019. [16] S. B. Fuller. Four Wings: An Insect-Sized Aerial Robot with Steering Ability and Payload Capacity for Autonomy. IEEE Robotics and Automation Letters, 4(2):570{577, 2019. [17] R. J. Wood. The First Takeo of a Biologically Inspired At-Scale Robotic Insect. IEEE Transactions on Robotics, 24(2):341{347, 2008. [18] X. Yang, Y. Chen, L Chang, A. A. Calder on, and N. O. P erez-Arancibia. Bee+: A 95-mg Four-Winged Insect-Scale Flying Robot Driven by Twinned Unimorph Actuators. IEEE Robotics and Automation Letters, 4:4270{4277, 2019. [19] Z. Tu, F. Fei, J. Zhang, and X. Deng. Acting Is Seeing: Navigating Tight Space Using Flapping Wings. arXiv preprint arXiv:1902.08688, 2019. [20] S. P. Sane and M. H. Dickinson. The Aerodynamic Eects of Wing Rotation and a Revised Quasi-Steady Model of Flapping Flight. Journal of Experimental Biology, 205(8):1087{1096, 2002. [21] W. B. Dickson and M. H. Dickinson. The Eect of Advance Ratio on the Aero- dynamics of Revolving Wings. Journal of Experimental Biology, 207(24):4269{ 4281, 2004. [22] D. Ishihara, T. Horie, and T. Niho. An Experimental and Three-Dimensional Computational Study on the Aerodynamic Contribution to the Passive Pitching Motion of Flapping Wings in Hovering Flies. Bioinspiration & Biomimetics, 9(4):046009, 2014. 127 [23] L. Chang and N. O. P erez-Arancibia. The Dynamics of Passive Wing-Pitching in Hovering Flight of Flapping Micro Air Vehicles Using Three-Dimensional Aerodynamic Simulations. In AIAA Atmospheric Flight Mechanics Conference, page 0013, San Diego, CA, 2016. [24] A. J. Bergou, S. Xu, and Z. J. Wang. Passive Wing Pitch Reversal in Insect Flight. Journal of Fluid Mechanics, 591:321{337, 2007. [25] Y. Chen, N. Gravish, A. L. Desbiens, R. Malka, and R. J. Wood. Experimental and Computational Studies of the Aerodynamic Performance of a Flapping and Passively Rotating Insect Wing. Journal of Fluid Mechanics, 791:1{33, 2016. [26] G. K. Taylor and A. L. R. Thomas. Animal Flight Dynamics II. Longitudinal Stability in Flapping Flight. Journal of Theoretical Biology, 214(3):351{370, 2002. [27] M. Sun, J. Wang, and Y. Xiong. Dynamic Flight Stability of Hovering Insects. Acta Mechanica Sinica, 23(3):231{246, 2007. [28] W. B. Dickson, A. D. Straw, and M. H. Dickinson. Integrative Model of Drosophila Flight. AIAA Journal, 46(9):2150{2164, 2008. [29] M. A. Bolender. Rigid Multi-Body Equations-of-Motion for Flapping Wing MAVs Using Kane's Equations. In AIAA Guidance, Navigation, and Control Conference, page 6158, Chicago, IL, 2009. [30] D. B. Doman, M. W. Oppenheimer, and D. O. Sigthorsson. Wingbeat Shape Modulation for Flapping-Wing Micro-Air-Vehicle Control During Hover. Jour- nal of Guidance, Control, and Dynamics, 33(3):724{739, 2010. [31] C. T. Orlowski and A. R. Girard. Dynamics, stability, and control analyses of apping wing micro-air vehicles. Progress in Aerospace Sciences, 51(Supplement C):18 { 30, 2012. [32] M. Sun. Insect Flight Dynamics: Stability and Control. Reviews of Modern Physics, 86(2):615, 2014. [33] J. K. Kim and J. Han. A Multibody Approach for 6-DOF Flight Dynamics and Stability Analysis of the Hawkmoth Manduca sexta. Bioinspiration & Biomimetics, 9(1):016011, 2014. [34] T. S. Clawson, S. B. Fuller, R. J. Wood, and S. Ferrari. A Blade Element Approach to Modeling Aerodynamic Flight of an Insect-Scale Robot. In 2017 American Control Conference (ACC), pages 2843{2849. IEEE, 2017. 128 [35] Z. J. Wang, J. M. Birch, and M. H. Dickinson. Unsteady Forces and Flows in Low Reynolds Number Hovering Flight: Two-Dimensional Computations vs Robotic Wing Experiments. Journal of Experimental Biology, 207(3):449{460, 2004. [36] J. H. Wu and M. Sun. Unsteady Aerodynamic Forces of a Flapping Wing. Journal of Experimental Biology, 207(7):1137{1150, 2004. [37] N. O. P erez-Arancibia, K. Y. Ma, K. C. Galloway, J. D. Greenberg, and R. J. Wood. First Controlled Vertical Flight of a Biologically Inspired Microrobot. Bioinspiration & Biomimetics, 6(3):036009, 2011. [38] B. M. Finio, N. O. P erez-Arancibia, and R. J. Wood. System Identication and Linear Time-Invariant Modeling of an Insect-Sized Flapping-Wing Micro Air Vehicle. In Intelligent Robots and Systems (IROS), 2011 IEEE/RSJ Interna- tional Conference on, pages 1107{1114. IEEE, 2011. [39] J. Wittenburg. Dynamics of Multibody Systems. Springer, 2007. [40] H. Baruh. Analytical Dynamics. WCB/McGraw-Hill Boston, 1999. [41] B. Liang and M. Sun. Nonlinear Flight Dynamics and Stability of Hovering Model Insects. Journal of The Royal Society Interface, 10(85):20130269, 2013. [42] S. P. Sane and M. H. Dickinson. The Control of Flight Force by a Flap- ping Wing: Lift and Drag Production. Journal of Experimental Biology, 204(15):2607{2626, 2001. [43] Q. T. Truong, Q. V. Nguyen, V. T. Truong, H. C. Park, D. Y. Byun, and N. S. Goo. A Modied Blade Element Theory for Estimation of Forces Generated by a Beetle-Mimicking Flapping Wing System. Bioinspiration & Biomimetics, 6(3):036008, 2011. [44] T. Nakata, H. Liu, and R. J. Bomphrey. A CFD-Informed Quasi-Steady Model of Flapping-Wing Aerodynamics. Journal of Fluid Mechanics, 783:323{343, 2015. [45] Y.J. Lee, K.B. Lua, T.T. Lim, and K.S. Yeo. A Quasi-Steady Aerodynamic Model for Flapping Flight with Improved Adaptability. Bioinspiration & Biomimetics, 11(3):036005, 2016. [46] C. P. Ellington. The Aerodynamics of Hovering Insect Flight. I. The Quasi- Steady Analysis. Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences, 305(1122):1{15, 1984. 129 [47] J. S. Han, J. K. Kim, J. W. Chang, and J. H. Han. An Improved Quasi-Steady Aerodynamic Model for Insect Wings that Considers Movement of the Center of Pressure. Bioinspiration & Biomimetics, 10(4):046014, 2015. [48] Edward C Polhamus. A concept of the vortex lift of sharp-edge delta wings based on a leading-edge-suction analogy. NASA Technical Note, 1966. [49] X. Deng, L. Schenato, W. C. Wu, and S. S. Sastry. Flapping Flight for Biomimetic Robotic Insects: Part I-System Modeling. IEEE Transactions on Robotics, 22(4):776{788, 2006. [50] A. J. Bergou, L. Ristroph, J. Guckenheimer, I. Cohen, and Z. J. Wang. Fruit Flies Modulate Passive Wing Pitching to Generate In-Flight Turns. Physical Review Letters, 104(14):148101, 2010. [51] S. E. Spagnolie, L. Moret, M. J. Shelley, and J. Zhang. Surprising Behaviors in Flapping Locomotion with Passive Pitching. Physics of Fluids, 22(4):041903, 2010. [52] H. Dai, H. Luo, and J. F. Doyle. Dynamic Pitching of an Elastic Rectangular Wing in Hovering Motion. Journal Fluid Mechanics, 693:473, 2012. [53] C. Kang and W. Shyy. Modeling of Instantaneous Passive Pitch of Flexible Flapping Wings. In 43rd AIAA Fluid Dynanmics Conference, page 2469, 2013. [54] F. Fei, J. A. Roll, and X. Deng. Design Principle of Wing Rotational Hinge Stiness. In 2015 IEEE International Conference on Robotics and Automation (ICRA), pages 1049{1054. IEEE, 2015. [55] H. Baruh. Another Look at the Describing Equations of Dynamics. The Journal of the Chinese Society of Mechanical Engineers, 21(1):15{23, 2000. [56] P. J. Dhrymes. Mathematics for Econometrics. Springer, 2013. [57] F. E. Udwadia and R. E. Kalaba. Analytical Dynamics: a New Approach. Cambridge University Press, 2007. [58] D. T. Greenwood. Advanced Dynamics. Cambridge University Press, 2006. [59] T. R. Kane and D. A. Levinson. Dynamics: Theory and Applications. McGraw Hill, 1985. [60] G. R. Kneller and K. Hinsen. Generalized Euler Equations for Linked Rigid Bodies. Physical Review E, 50(2):1559, 1994. [61] F. E. Udwadia. Introdcution to Quaternions and Their Application to Rigid Body Dynamics. USC, 2007. 130 [62] F. E. Udwadia and A. D. Schutte. An alternative derivation of the quaternion equations of motion for rigid-body rotational dynamics. Journal of Applied Mechanics, 77(4):044505, 2010. [63] Y. Xiong and M. Sun. Dynamic Flight Stability of a Bumblebee in Forward Flight. Acta Mechanica Sinica, 24(1):25{36, 2008. [64] H. E Taha, M. R. Hajj, and A. H. Nayfeh. Flight Dynamics and Control of Flapping-Wing MAVs: a Review. Nonlinear Dynamics, 70(2):907{939, 2012. [65] P. Chirarattananon, Y. Chen, E. F. Helbling, K. Y. Ma, R. Cheng, and R. J. Wood. Dynamics and Flight Control of a Flapping-Wing Robotic Insect in the Presence of Wind Gusts. Interface Focus, 7(1):20160080, 2017. [66] L. Chang and N. O. P erez-Arancibia. Time-Averaged Dynamic Modeling of A Flapping-Wing Micro Air Vehicle with Passive Rotation Mechanisms. In 2018 Atmospheric Flight Mechanics Conference, page 2830, 2018. [67] N. O. P erez-Arancibia, P. J. Duhamel, K. Y. Ma, and R. J. Wood. Model-Free Control of a Hovering Flapping-Wing Microrobot. Journal of Intelligent and Robotic Systems, 77(1):95{111, 2015. [68] K. Y. Ma, P. Chirarattananon, S. B. Fuller, and R. J. Wood. Controlled Flight of a Biologically Inspired, Insect-Scale Robot. Science, 340(6132):603{607, 2013. [69] H. Nagai, K. Isogai, T. Fujimoto, and T. Hayase. Experimental and Numerical Study of Forward Flight Aerodynamics of Insect Flapping Wing. AIAA Journal, 47(3):730{742, 2009. [70] J. S. Han, J. W. Chang, and J. H. Han. An Aerodynamic Model for Insect Flap- ping Wings in Forward Flight. Bioinspiration & Biomimetics, 12(3):036004, 2017. [71] M. S. Pantazopoulos. Vortex-Induced Vibration Parameters: Critical Review. In International Conference on Oshore Mechanics and Arctic Engineering (OMAE), volume 1, pages 199{255, 1994. [72] S. B. Fuller, Z. E. Teoh, P. Chirarattananon, N. O. P erez-Arancibia, J. Green- berg, and R. J. Wood. Stabilizing Air Dampers for Hovering Aerial Robotics: Design, Insect-Scale Flight Tests, and Scaling. Autonomous Robots, 41(8):1555{ 1573, 2017. [73] N. Gravish and R. J. Wood. Anomalous Yaw Torque Generation from Pas- sively Pitching Wings. In 2016 IEEE International Conference on Robotics and Automation (ICRA), pages 3282{3287. IEEE, 2016. 131 [74] Y. Wang, X. Yang, Y. Chen, D. K. Wainwright, C. P. Kenaley, Z. Gong, Z. Liu, H. Liu, J. Guan, T. Wang, J. C. Weaver, R. J. Wood, and Li Wen. A biorobotic adhesive disc for underwater hitchhiking inspired by the remora suckersh. Sci- ence Robotics, 2(10):eaan8072, 2017. [75] R. J. Wood. Design, fabrication, and analysis of a 3dof, 3cm apping-wing mav. In 2007 IEEE/RSJ international conference on intelligent robots and systems, pages 1576{1581. IEEE, 2007. [76] N. T. Jaeris, M. Lok, N. Winey, G.Y. Wei, and R. J. Wood. Multilayer lami- nated piezoelectric bending actuators: Design and manufacturing for optimum power density and eciency. Smart Materials and Structures, 25(5):055033, 2016. [77] NASA. User's Manual for OVERFLOW 2.2. https://overflow.larc.nasa. gov/home/users-manual-for-overflow-2-2/, 2019. [78] William D. H. A Primer for Wring PDE Solvers with Overture. https:// overtureframework.org/, 2011. [79] ANL. NEK5000. https://nek5000.mcs.anl.gov/, 2011. [80] ESI CFD. OpenFOAM User Guide. https://openfoam.com/, 2019. [81] T. D. Economon, F. Palacios, S. R. Copeland, T. W. Lukaczyk, and J. J. Alonso. SU2: An Open-Source Suite for Multiphysics Simulation and Design. AIAA Journal, 54(3):828{846, 2015. [82] H. Hadzic. Development and Application of Finite Volume Method for the Com- putation Of Flows around Moving Bodies on Unstructured, Overlapping Grids. PhD thesis, Technische Universit at Hamburg, 2006. [83] ANSYS, Inc. ANSYS Fluent Theory Guide, 2013. [84] R. Mittal and G. Iaccarino. Immersed Boundary Methods. Annu. Rev. Fluid Mech., 37:239{261, 2005. [85] C. S. Peskin. The Immersed Boundary Method. Acta Numerica, 11:479{517, 2002. [86] H. S. Udaykumar, R. Mittal, P. Rampunggoon, and A. Khanna. A Sharp Interface Cartesian Grid Method for Simulating Flows with Complex Moving Boundaries. Journal of Computational Physics, 174(1):345{380, 2001. [87] J. Donea, A. Huerta, J.P. Ponthot, and A. Rodr guez-Ferran. Arbitrary Lagrangian{Eulerian Methods. Encyclopedia of Computational Mechanics Sec- ond Edition, pages 1{23, 2017. 132 [88] CD-adapco. STAR-CCM+ Version 12.02 User Guide, 2017. [89] W. M. Chan. Overset Grid Technology Development at Nasa Ames Research Center. Computers & Fluids, 38(3):496{503, 2009. [90] W. D. Henshaw. Adaptive Mesh and Overlapping Grid Methods. Encyclopedia of Aerospace Engineering, 2010. [91] H. S. Tang, S. C. Jones, and F. Sotiropoulos. An Overset-Grid Method for 3D Unsteady Incompressible Flows. Journal of Computational Physics, 191(2):567{ 600, 2003. [92] J. Cai, H. M. Tsai, and F. Liu. A Parallel Viscous Flow Solver on Multi-Block Overset Grids. Computers & Fluids, 35(10):1290{1301, 2006. [93] J. McNaughton, I. Afgan, D. D. Apsley, S. Rolfo, T. Stallard, and P. K. Stansby. A Simple Sliding-Mesh Interface Procedure and Its Application to the Cfd Simu- lation of a Tidal-Stream Turbine. International Journal for Numerical Methods in Fluids, 74(4):250{269, 2014. [94] R. Lohner. Robust, Vectorized Search Algorithms for Interpolation on Unstruc- tured Grids. Journal of Computational Physics, 118(2):380{387, 1995. [95] J. H. Ferziger and M. Peric. Computational Methods for Fluid Dynamics. Springer Science & Business Media, 2012. [96] H. K. Versteeg and W. Malalasekera. An Introduction to Computational Fluid Dynamics: The Finite Volume Method. Pearson education, 2007. [97] S. R. Mathur and J. Y. Murthy. A Pressure-based Method for Unstructured Meshes. Numerical Heat Transfer, 31(2):195{215, 1997. [98] S. Patankar. Numerical Heat Transfer and Fluid Flow. CRC press, 2018. [99] I. Demird zi c and M. Peri c. Space Conservation Law in Finite Volume Calcu- lations of Fluid Flow. International Journal for Numerical Methods in Fluids, 8(9):1037{1050, 1988. [100] Manchester CFD group. Finite Volume Methods for Unstructured Grids, 2011. [101] B. Yu, W.Q. Tao, J.J. Wei, Y. Kawaguchi, T. Tagawa, and H. Ozoe. Discussion on Momentum Interpolation Method for Collocated Grids of Incompressible Flow. Numerical Heat Transfer: Part B: Fundamentals, 42(2):141{166, 2002. [102] D. K. Kolmogorov, W. Z. Shen, N. N. Srensen, and J. N. Srensen. Fully Con- sistent SIMPLE-like algorithms on collocated Grids. Numerical Heat Transfer, Part B: Fundamentals, 67(2):101{123, 2015. 133 [103] A. Pascau. Cell Face Velocity Alternatives in A Structured Colocated Grid for the Unsteady Navier{Stokes Equations. International Journal for Numerical Methods in Fluids, 65(7):812{833, 2011. [104] W. Z. Shen, J. A. Michelsen, and J. N. Sorensen. Improved Rhie-Chow Inter- polation for Unsteady Flow Computations. AIAA Journal, 39(12):2406{2409, 2001. 134 Appendices Appendix A. Matrices and derivatives The rotation matrices for the body and wing i are R b = 2 4 c b c b c b s b s b s b s b c b c b s b s b s b s b +c b c b s b c b c b s b c b +s b s b c b s b s b s b c b c b c b 3 5 ; (A.1a) R i = 2 4 c i c i c i s i c i +s i s i c i s i s i s i c i s i c i c i c i s i s i c i s i s i c i c i s i s i s i s i +c i c i 3 5 ; (A.1b) where c = cos(), s = sin(). In this study, the deviation angle of the wings is assumed to be zero during ight, i.e., i = 0. The proof of the equality _ R T = R T ~ ! is available in [57]. For a column vector = [ 1 ; 2 ; 3 ] T , the skew symmetric matrix is dened as ~ = 2 4 0 3 2 3 0 1 2 1 0 3 5 (A.2) The partial derivatives of angular and translational velocities with respect to _ q are @! b @ _ q = [0; T b ; 0; 0]; @v b @ _ q = [I; 0; 0; 0] (A.3a) @! 1 @ _ q = [0; R 1 T b ; T 1 ; 0]; @! 2 @ _ q = [0; R 2 T b ; 0; T 2 ]; (A.3b) @v 1 @ _ q = I;R T b ~ r r;1 T b R T b R T 1 ~ r c;1 R 1 T b ;R T b R T 1 ~ r c;1 T 1 ; 0 ; (A.3c) @v 2 @ _ q = I;R T b ~ r r;2 T b R T b R T 2 ~ r c;2 R 2 T b ; 0;R T b R T 2 ~ r c;2 T 2 ; (A.3d) @v o r;i @ _ q = I;R T b ~ r r;i T b ; 0; 0 (A.3e) 135 The time rates of the partial derivatives of the angular velocities are d dt @! b @ _ q = h 0; _ T b ; 0; 0 i ; (A.4a) d dt @! 1 @ _ q = h 0; _ R 1 T b + R 1 _ T b ; _ T 1 ; 0 i ; (A.4b) d dt @! 2 @ _ q = h 0; _ R 2 T b + R 2 _ T b ; 0; _ T 2 i : (A.4c) The partial derivatives of angular velocities with respect to q are @! b @q = 0; @T b @ b _ b ; 0; 0 (A.5a) @! 1 @q = 0; R 1 @T b @ b _ b ; @R 1 @ 1 ! b ( 1 ) + @T 1 @ 1 _ 1 ; 0 ; (A.5b) @! 2 @q = 0; R 2 @T b @ b _ b ; 0; @R 2 @ 2 ! b ( 2 ) + @T 2 @ 2 _ 2 ; (A.5c) The partial derivative of the matrix T with respect to the coordinate vector b is dened as @T b @ b = @T b @ b ; @T b @ b ; @T b @ b : (A.6) For simplicity, this denition is dierent from the one in [56]. The same denition is applied to @R i @ i and @T i @ i . The associated matrix operator () is dened as () = 2 6 6 6 4 0 0 0 0 . . . . . . . . . . . . 0 0 3 7 7 7 5 s ; (A.7) where is as1 column vector and is a 1 column vector. This operator denes a block diagonal matrix with and. If the vectors and have the same dimension, the parenthesis part is ignored in the notation. The time derivatives of the angular and translational velocities are d! b dt = _ T b _ b + T b b ; dv b dt = _ v b ; (A.8a) d! i dt = _ R i T b + R i _ T b _ b + R i T b b + _ T i _ i + T i i ; (A.8b) 136 dv i dt = _ v b _ R T b ~ r r;i ! b R T b ~ r c;i _ ! b _ R T b R T i + R T b _ R T i ~ r c;i ! i R T b R T i ~ r c;i _ ! i : (A.8c) The mass matrix is formed by four blocks corresponding to the generalized veloc- ities _ v b , b , 1 , 2 terms M = [M 1 ; M 2 ; M 3 ; M 4 ]; (A.9a) where the partition matrices are M 1 = @v b @ _ q T m b + 2 X i=1 @v i @ _ q T m i ; (A.9b) M 2 = @! b @ _ q T J b T b + 2 X i=1 @! i @ _ q T J i R i T b + 2 X i=1 @v i @ _ q T m i R T b ~ r r;i T b R T b R T i ~ r c;i R i T b ; (A.9c) M 3 = @v 1 @ _ q T m 1 R T b R T 1 ~ r c;1 T 1 + @! 1 @ _ q T J 1 T 1 ; (A.9d) M 4 = @v 2 @ _ q T m 2 R T b R T 2 ~ r c;2 T 2 + @! 2 @ _ q T J 2 T 2 : (A.9e) Similarly, the structural damping matrix is dened as C = [C 1 ; C 2 ; C 3 ; C 4 ]; (A.10a) where the partition matrices are C 1 = 0; (A.10b) C 2 = @! b @ _ q T J b _ T b + d dt @! b @ _ q @! b @q T J b T b + 2 X i=1 @v i @ _ q T m i n _ R T b ~ r r;i T b R T b ~ r r;i _ T b _ R T b R T i + R T b _ R T i ~ r c;i R i T b R T b R T i ~ r c;i _ R i T b + R i _ T b o + 2 X i=1 @! i @ _ q T J i _ R i T b + R i _ T b + d dt @! i @ _ q @! i @q T J i R i T b ; (A.10c) 137 C 3 = @v 1 @ _ q T m 1 n _ R T b R T 1 + R T b _ R T 1 ~ r c;1 T 1 R T b R T 1 ~ r c;1 _ T 1 o + @! 1 @ _ q T J 1 _ T 1 + d dt @! 1 @ _ q @! 1 @q T J 1 T 1 ; (A.10d) C 4 = @v 2 @ _ q T m 2 n _ R T b R T 2 + R T b _ R T 2 ~ r c;2 T 2 R T b R T 2 ~ r c;2 _ T 2 o + @! 2 @ _ q T J 2 _ T 2 + d dt @! 2 @ _ q @! 2 @q T J 2 T 2 : (A.10e) The coordinate vector of the potential energy derivative is @V @q T = 2 6 6 6 6 6 6 6 6 6 6 4 (m b +m 1 +m 2 )gu 2 X i=1 m i g r r;i + R T i r c;i T @R T b @ b T u m 1 g r T c;1 ( 1 ) @R T 1 @ 1 T R b u + [0;k h;1 1 ] T m 2 g r T c;2 ( 2 ) @R T 2 @ 2 T R b u + [0;k h;2 2 ] T 3 7 7 7 7 7 7 7 7 7 7 5 : (A.11) The matrices for the reduced order model are M r = 2 4 M (1:6;1:6) M (1:6;8) M (1:6;10) M (8;1:6) M (8;8) M (8;10) M (10;1:6) M (10;8) M (10;10) 3 5 ; C r = 2 4 C (1:6;1:6) C (1:6;8) C (1:6;10) C (8;1:6) C (8;8) C (8;10) C (10;1:6) C (10;8) C (10;10) : 3 5 ; D r = 2 4 D (1:6) D (8) D (10) 3 5 ;A i = 2 4 M (1:6;j) M (8;j) M (8;j) 3 5 ;B i = 2 4 C (1:6; j) C (8;j) C (10;j) 3 5 ; (A.12) where j = 7 + 2(i 1), i=1,2. 138 Appendix B. Parameters of the wing and body Geometric and inertial parameters for the body and wings are required for both the CFD and complete dynamics simulations. The characteristic dimensions of the wing for MAV I are shown in Fig. A.1(a). The leading edge and wing root are the geometric reference locations for the wing. In Fig. A.1(b), the top of the head and longitudinal symmetric plane are the geometric reference locations for the body. The centers of mass of the wing and body are chosen to be their respective inertial reference points. The shapes of the wing and body for MAV II are plotted in Figs. A.1(c)(d). The values of parameters for the two robotic MAVs are listed in Table A.1-A.4. (a) (b) (c) (d) x c,i y c,i x r,i y r,i l s l h c max l b l c Fig. A.1. Geometry and diagram of the wings and body for the two MAVs. (a)-(b) MAV I. (c)-(d) MAV II. Table A.1: Geometric parameters for the wings. MAV l s c max c l h x c,2 y c,1 x c,2 y c,2 I 15:3 4.0 2.9 1.3 0.2841 6.902 0.2841 6.902 II 15.4 5.7 3.4 6.4 0.837 6.041 0.837 6.041 Notes: Units for length, mass, moment of inertia and hinge stiness are mm, mg, mg mm 2 , Nm, respectively. This convention is applied to Tables A.1A.4. Table A.2: Inertial parameters for the wings. MAV m 1 J 1 (1;1) J 1 (1;2) J 1 (1;3) J 1 (2;2) J 1 (2;3) J 1 (3;3) t J 1 (2;2) t J 1 (2;3) I 1:192 16.989 -0.261 0.0 1.268 0.0 18.254 1.365 -2.599 II 1.09 17.02 -0.31 0.0 2.16 0 19.13 2.882 -5.84 139 Table A.3: Geometric parameters for the bodies. MAV l b l c x r,1 y r,1 z r,1 x r,2 y r,2 z r,2 I 27:3 14.5 6.554 2.5 0.0 6.554 2.5 0.0 II 14.0 9 8.15 1.235 -0.1 8.15 1.235 -0.1 Table A.4: Inertial parameters for the bodies. MAV m b J b (1;1) J b (1;2) J b (1;3) J b (2;2) J b (2;3) J b (3;3) I 60:0 215.553 0.0 0.0 4443.85 0.0 4443.85 II 69.0 67.0 0.0 21.3 1272.3 0.0 1286.3 140 Appendix C. File locking for parallel implementa- tion of User Code 1 2 void EqnDyn( double Time , double TimeStep , i n t I t e r ) f 3 4 i n t i , j ; 5 double t ; 6 double tn ; 7 double dt ; 8 FILE fp ; 9 i n t f i l e d e s c ; 10 11 char d i r l c k [ 2 5 5 ] = " lck . f i l e " ; 12 char diroutput [ 2 5 5 ] = " output . dat " ; 13 char dirinput [ 2 5 5 ] = " input . csv " ; 14 15 t = Time ; 16 dt = TimeStep ; 17 18 while (1) f 19 f i l e d e s c = open ( d i r l c k , O RDWR j O CREAT j O EXCL, 0444); 20 i f ( f i l e d e s c != 1) f 21 fp = fopen ( diroutput , " r " ) ; 22 f s c a n f ( fp , "%l f " , &tn ) ; 23 f c l o s e ( fp ) ; 24 25 i f ( fabs ( t tn ) > 0.5 dt ) f 26 fp = fopen ( dirinput , "a+" ) ; 27 f p r i n t f ( fp , "%.6 l f ,%.6 l f , %dnn" , t , dt , I t e r ) ; 28 f c l o s e ( fp ) ; 29 30 system ( "matlab nodesktop nosplash r n" run 31 SolveEqns ; e x i t ;n" " ) ; 32 33 usleep ( 1 0 0 ) ; // writing of output . dat in SolveEqns .m 34 35 g 36 37 ( void ) c l o s e ( f i l e d e s c ) ; 38 ( void ) unlink ( d i r l c k ) ; 39 40 break ; 41 g 42 43 g 44 45 g 141

Abstract (if available)

Abstract Birds and insects inspired us to develop flying machines, such as airplanes and helicopters. However, the development of artificial flapping-wing flyers falls far behind that of other aircraft, which is possibly due to the complexities of the associated unsteady aerodynamics, bidirectional actuation, stability and control, transmission design, etc. To aid the practical design of insect-sized robotic micro air vehicles (MAVs), we aimed to develop a set of tools to analyze the unsteady aerodynamics and flight dynamics of flapping-wing flyers. ❧ The robotic design of MAVs with wings connected to the body through flexure hinges is able to reduce the number of actuation degrees of freedom and the complexity of transmission design. In this design, each wing is partially actuated and the pitching motion of the wing is dynamically determined under the interactions among inertial, hinge elastic, and aerodynamic forces. Both the wing motion itself and the dynamics of the entire MAV define fluid-structure interaction (FSI) procedures. ❧ For estimating the aerodynamic forces on a single wing in passive pitching, the dynamics model of one wing is first derived and simulated by integrating the equations of motion with a proprietary computational fluid dynamics (CFD) solver. To achieve physical insights and quick force estimation, revised empirical formulas are proposed to extend the adaptability and improve the accuracy of the existing quasi-steady (QS) models. Preliminary models for the pitching angle limit applied by hinge geometry are proposed and simulated. The validity of the weakly decoupled FSI assumption is discussed with the simulation of a single wing immersed in water. ❧ The full dynamics model of the entire MAV is formulated using our proposed alternative form of the Lagrange's equations. In this alternative form, repeated inertial terms are removed to enable the manual derivation of compact equations of motion for the MAV to achieve efficient numerical computations. The integrated dynamics simulation with the proprietary CFD solver provides us with high-fidelity numerical predictions for the dynamic behavior of the flapping-wing MAVs in open-loop or closed-loop simulations. The computationally efficient full dynamics model is formulated with QS formulas. Simulation results are compared with the results by the integrated simulation with CFD solver. The full dynamics model with QS formulas is experimentally validated using both the underwater and in-air jump tests for quantitative and qualitative purposes, respectively. Highly satisfactory matching between the simulated and experimental results is achieved. The full dynamics model is proved to be effective in the dynamic behavior prediction for the robotic MAVs. ❧ To compensate for the incapability of the proprietary CFD solver in flexibility and computing efficiency, a unified CFD solver based on the overset grid method is proposed and implemented for moving boundary problems at low Reynolds numbers. To achieve efficient and accurate simulations with the segregated overset grid solver, a fast search-based hole-cutting algorithm, a hybrid automatic multigrid generation process, and several consistency-improved discretization schemes are proposed. Finally, testing examples are provided to demonstrate the capability of the in-house solver on unsteady and moving boundary problems.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Wing selection and design for flapping flight at small scale

PDF

Passive rolling and flapping dynamics of a heaving Λ ﬂyer

PDF

Novel soft and micro transducers for biologically-inspired robots

PDF

Development and applications of a body-force propulsor model for high-fidelity CFD

PDF

Development of biologically-inspired sub-gram insect-scale autonomous robots

PDF

Study of belly-flaps to enhance lift and pitching moment coefficient of a blended-wing-body airplane in landing and takeoff configuration

PDF

The projection immersed boundary method for compressible flow and its application in fluid-structure interaction simulations of parachute systems

PDF

Mechanical and flow interactions facilitate cooperative transport and collective locomotion in animal groups

PDF

BeeCAS: a collision avoidance system for micro-aerial vehicles

PDF

Numerical study of shock-wave/turbulent boundary layer interactions on flexible and rigid panels with wall-modeled large-eddy simulations

PDF

Simulation and machine learning at exascale

PDF

Development and control of biologically-inspired robots driven by artificial muscles

PDF

Modeling and simulation of complex recovery processes

PDF

Numerical modeling of separated flows at moderate Reynolds numbers appropriate for turbine blades and unmanned aero vehicles

PDF

Evaluating sensing and control in underwater animal behaviors

PDF

The effect of lattice structure and porosity on thermal conductivity of additively-manufactured porous materials

PDF

Plant substructuring and real-time simulation using model reduction

PDF

Computationally efficient design of optimal strategies for passive and semiactive damping devices in smart structures

PDF

Modeling of bicontinuous nanoporous materials and investigation of their mechanical and hydraulic properties

PDF

High throughput computational framework for synthesis and accelerated discovery of dielectric polymer materials using polarizable reactive molecular dynamics and graph neural networks

Asset Metadata

Creator Chang, Longlong (author)

Core Title Dynamic modeling and simulation of flapping-wing micro air vehicles

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Aerospace Engineering

Publication Date 12/16/2019

Defense Date 09/20/2019

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag computational fluid dynamics,dynamic modeling,flapping wing,fluid-structure interaction,micro air vehicles,OAI-PMH Harvest

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Perez-Arancibia, Nestor (committee chair), Domaradzki, Julian (committee member), Luhar, Mitul (committee member), Nakano, Aiichiro (committee member)

Creator Email longlonc@usc.edu,sdcll221@gmail.com

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c89-254493

Unique identifier UC11673437

Identifier etd-ChangLongl-8046.pdf (filename),usctheses-c89-254493 (legacy record id)

Legacy Identifier etd-ChangLongl-8046.pdf

Dmrecord 254493

Document Type Dissertation

Rights Chang, Longlong

Type texts

Source 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 a...

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