<%BANNER%>

Flapping and Fixed Wing Aerodynamics of Low Reynolds Number Flight Vehicles


PAGE 1

FLAPPING AND FIXED WING AERODYNAMICS OF LOW REYNOLDS NUMBE R FLIGHT VEHICLES By DRAGOS VIIERU A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLOR IDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2006

PAGE 2

Copyright 2004 by Dragos Viieru

PAGE 3

ACKNOWLEDGMENTS I would like to express my deepest appreciation to my advisor, Dr. Wei Shyy, for his guidance and support. I am also grateful for his patience and constant encouragement throughout my research activity. I owe special thanks to Dr. Peter Ifju and his micro air vehicles research group for sharing their knowledge and experience. I would also like to thank Dr. Renwei Mei, Dr. Corin Segal and Dr. William Hager for serving as members of my supervisory committee. I acknowledge the close collaboration with Dr. Roberto Albertani regarding experimental aspects of micro air vehicle aerodynamics and Kyuho Lee who kindly provided me with his wing shape designs that made part of the work in chapter 3 possible. I would like to thank all my colleagues of the computational thermo-fluids group for their friendship and support, in particular, Dr. Siddharth Thakur and Dr. Ramji Kamakoti for generously sharing their time to answer my queries. I would like to thank Dr. Yongsheng Lian for his collaboration during my first year at University of Florida. My thanks also go to Dr. Jian Tang, who shared his broad experience on flapping flight aerodynamics. Finally, I would like to thank to my parents for their understanding, patience and encouragement when it was most needed. iii

PAGE 4

TABLE OF CONTENTS page ACKNOWLEDGMENTS.................................................................................................iii LIST OF TABLES.............................................................................................................vi LIST OF FIGURES..........................................................................................................vii ABSTRACT.........................................................................................................................x CHAPTER 1. INTRODUCTION TO FLAPPING FLIGHT..............................................................1 1.1 Introduction.............................................................................................................1 1.2 Flapping Wing in Nature........................................................................................2 1.2.1 Wing Motions...............................................................................................3 1.2.2 Flight Modes in Birds...................................................................................4 1.2.3 Beating Frequency and Scaling....................................................................5 1.3 Insect Flight..........................................................................................................10 1.3.1 General Characteristics...............................................................................10 1.3.2 Experimental Efforts in Studying the Insect Flight Mechanisms...............11 1.3.3 Computational Efforts in Studying the Insect Flight Mechanisms............13 1.4 Unsteady mechanisms identified to enhance lift in insects..................................18 1.4.1 Analytical and quasi-steady models of insect flight...................................18 1.4.2 Wagner Effect.............................................................................................20 1.4.3 Leading-edge Vortex and Delayed Stall.....................................................21 1.4.4 Rotational Forces........................................................................................24 1.4.5 Wing-wake Interaction...............................................................................25 1.4.6 Clap and Fling Mechanism.........................................................................27 1.5 Objective and Outline...........................................................................................28 2. COMPUTATIONAL METHODS FOR COUPLED FLUID-STRUCTURE INTERACTIONS.......................................................................................................30 2.1 Overview...............................................................................................................30 2.2 Governing Equations for Fluid Flows..................................................................31 2.3 Fluid Flow Solvers................................................................................................33 2.4 Geometric Conservation Law...............................................................................34 iv

PAGE 5

2.5 Moving Grid Technique.......................................................................................35 2.5.1 Master-slave Concept.................................................................................35 2.5.2 Grid Re-arrangement Algorithm with Weighting Functions......................39 2.5.3 Grid Re-arrangement Algorithm Based on Spring Analogy Concept........41 2.6 Membrane Structural Solver.................................................................................42 2.6.1 Mooney-Rivilin Model...............................................................................44 2.6.2 Solution of the Dynamic Equations............................................................45 2.7 Coupling between Fluid Solver and Structural Solver.........................................47 3. FIXED WING MICRO AIR VEHICLES..................................................................51 3.1 Overview...............................................................................................................51 3.2 Rigid Wing Micro Air Vehicles...........................................................................54 3.2.1 Wing Geometry..........................................................................................54 3.2.2 Computational Setup..................................................................................54 3.2.3 Tip Vortex Effect on the Rigid Wing Aerodynamics.................................56 3.3 Flexible Wing Micro Air Vehicles.......................................................................65 3.4 Summary and Conclusions for Fixed Wing Aerodynamics.................................70 4. AERODYNAMICS OF LOW REYNOLDS NUMBER HOVERING AIRFOILS...72 4.1 Validation of Viscous Force Computation...........................................................72 4.1.1 Steady Flow over a Flat Plate at Zero Angle of Attack..............................73 4.1.2 Flow over an Oscillating Flat Plate in a Quiescent Environment..............76 4.2 Kinematics of 2-D Airfoil in Hovering Flight......................................................79 4.2.1 Kinematics of Water Treading Mode.....................................................80 4.2.2 Kinematics of Normal Hovering Mode..................................................82 4.3 Water Treading Hovering Mode.......................................................................82 4.3.1 Grid and Time Increment Sensitivity Analysis..........................................84 4.3.2 Comparison with Similar Work..................................................................87 4.3.3 Lift Generation Mechanism for Water Treading Hovering Mode..........87 4.4 Normal Hovering Mode....................................................................................91 4.4.1 Grid Refinement Study...............................................................................91 4.4.2 Comparison with Similar Computational and Experimental Results.........93 4.5 Lift Generation Mechanisms in Normal and Water Treading Hovering Modes.....................................................................................................................96 4.6 Summary and Conclusion for the Aerodynamics of Hovering Airfoils.............104 5. SUMMARY AND FUTURE WORK......................................................................107 5.1 Present Progress of Research Work....................................................................107 5.2 Future Research Direction..................................................................................109 LIST OF REFERENCES.................................................................................................111 BIOGRAPHICAL SKETCH...........................................................................................124 v

PAGE 6

LIST OF TABLES Table page 1-1. Variety of beating frequency and Reynolds number for birds and insects...............7 1-2. Experimental robotic and mechanic devices for insect flapping wing study..........14 1-3. Principal numerical studies of insect flapping flight...............................................18 2-1. Moving boundaries techniques. Advantages and disadvantages.............................42 3-1. Aerodynamic performances for flexible and rigid wing.........................................69 4-1. Numerical and analytical friction coefficient values for flat plate..........................75 4-2. Grid and time increment details for the water treading hovering cases...............84 4-3. Grid and time increment details for the normal hovering cases..........................91 4-4. Kinematics parameters for the experimental and computational setup...................94 4-5. Kinematic parameters for water treading and normal hovering modes. The Reynolds number for both cases is 100...................................................................97 4-6. Parameters for water treading hovering mode employed for Reynolds number effect study............................................................................................................101 vi

PAGE 7

LIST OF FIGURES Figure page 1-1. Wing beating motions...............................................................................................3 1-2. Hovering flight..........................................................................................................5 1-3. Wing beat frequency versus the body mass for different insect species...................9 1-4. Cross-sectional corrugations of dragonfly wings....................................................12 1-5. Leading edge vortex on a hawkmoth wing.............................................................15 1-6. Flow around a thin airfoil........................................................................................22 1-7. Leading edge vortex (LEV) structures at different Reynolds numbers...................24 1-8. Momentum transfer due to wing-wake interaction.................................................26 1-9. Schematic of clap and fling mechanism..................................................................27 2-1. Illustration of the moving grid technique................................................................38 2-2. Master-slave algorithm performance for large rotations.........................................39 2-3. Current node neighbors position for grid re-arrangement algorithm......................40 2-4. Grid quality for the grid re-arrangement algorithm.................................................41 2-5. Generic algorithm for fluid-structure interaction....................................................49 3-1. Fixed and flapping wing MAVs..............................................................................53 3-2. MAVs wing shapes................................................................................................54 3-3. Computational domain for the MAV wing.............................................................55 3-4. Streamlines, vorticity and pressure on the 6 MAV wings upper surface as a function of the angle of attack.................................................................................57 3-5. Wing shape geometry..............................................................................................59 vii

PAGE 8

3-6. Vorticity magnitude contours behind the wing at 6 o angle of attack......................60 3-7. Pressure contours on the wings upper surface and streamlines slightly above the wing at 6 o angle of attack........................................................................................61 3-8. Pressure contours on the wings lower surface and streamlines slightly below the wing at 6 o angle of attack........................................................................................62 3-9. Pressure coefficient on the wing surface at x/c = 0.34 and 6 o angle of attack........63 3-10. Pressure coefficient on the wing surface at x/c = 0.53 and 6 o angle of attack........63 3-11. Span-wise lift distribution at 6 o angle of attack......................................................64 3-12. Span-wise drag distribution at 6 o angle of attack....................................................65 3-13. Surface grid for the 5 MAV wing..........................................................................66 3-14. History of the vertical displacements of two points, one near the wing tip and the other one between the two battens..........................................................................67 3-15. Vertical displacement contours at various times.....................................................68 3-16. Aerodynamic forces history for the 5 flexible wing MAV...................................68 4-1. Numerical and analytical velocity distribution in the boundary layer along different locations on a flat plate and different Reynolds numbers.......................................74 4-2. Numerical and analytical friction coefficient C f for a flat plate versus Reynolds number.....................................................................................................................76 4-3. Skin-friction and wall velocity for an oscillating plate...........................................77 4-4. Velocity distribution above an oscillating plate......................................................78 4-5. Hummingbird in hovering flight.............................................................................79 4-6. The schematics of two hovering styles....................................................................80 4-7. Schematics of the water treading hovering mode................................................81 4-8. Schematics of normal hovering mode.................................................................82 4-9. The coarse multi-block grid used in the current work. The grid has 90 points around the airfoil and 70 points in the radial direction...........................................83 4-10. Lift coefficient for three periods and different grid sizes and t=0.01....................85 4-11. Drag coefficient for three periods and different grid sizes and t=0.01..................85 viii

PAGE 9

4-12. Lift coefficient for three periods and different time increments.............................86 4-13. Drag coefficient for three periods and different time increments...........................86 4-14. Lift and Drag coefficients for one period................................................................88 4-15. Kinematics of water treading hovering mode and the aerodynamic forces for one complete stroke........................................................................................................89 4-16. Flow field snapshots of the water treading hovering mode.................................90 4-17. Lift coefficient of the normal hovering mode for three periods and different grid sizes and t=0.01.....................................................................................................92 4-18. Lift coefficient history for the normal hovering mode........................................95 4-19. Vorticity plot for normal hovering mode with h a /c=2.4 and Re=115..................96 4-20. One cycle force history for two hovering modes....................................................98 4-21. Vorticity contours for two hovering modes............................................................99 4-22. Lift coefficient for the water treading mode......................................................102 4-23. Vorticity contours for the water treading mode.................................................103 4-24. Pressure distribution on the airfoil surface for the water treading mode...........104 ix

PAGE 10

Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy FLAPPING AND FIXED WING AERODYNAMICS OF LOW REYNOLDS NUMBER FLIGHT VEHICLES By Dragos Viieru August 2006 Chair: Wei Shyy Major Department: Mechanical and Aerospace Engineering Lately, micro air vehicles (MAVs), with a maximum dimension of 15 cm and nominal flight speed around 10m/s, have attracted interest from scientific and engineering communities due to their potential to perform desirable flight missions and exhibit unconventional aerodynamics, control, and structural characteristics, compared to larger flight vehicles. Since MAVs operate at a Reynolds number of 10 5 or lower, the lift-to-drag ratio is noticeably lower than the larger manned flight vehicles. The light weight and low flight speed cause MAVs to be sensitive to wind gusts. The MAVs small overall dimensions result in low aspect ratio wings with strong wing tip vortices that further complicate the aerodynamics of such vehicles. In this work, two vehicle concepts are considered, namely, fixed wings with flexible structure aimed at passive shape control, and flapping wings aimed at enhancing aerodynamic performance using unsteady flow fields. A finite volume, pressure-based Navier-Stokes solver along with moving grid algorithms is employed to simulate the flow field. The coupled fluid-structural dynamics x

PAGE 11

of the flexible wing is treated using a hyperelastic finite element structural model, the above-mentioned fluid solver via the moving grid technique, and the geometric conservation law. Three dimensional aerodynamics around a low aspect ratio wing for both rigid and flexible structures and fluid-structure interactions for flexible structures have been investigated. In the Reynolds numbers range of 7x10 4 to 9x10 4 the flexible wing exhibits self-initiated vibrations even in steady free-stream, and is found to have a similar performance to the identical rigid wing for modest angles of attack. For flapping wings, efforts are made to improve our understanding of the unsteady fluid physics related to the lift generation mechanism at low Reynolds numbers (75 to 1,700). Alternative moving grid algorithms, capable of handling the large movements of the boundaries (characteristic of flapping wing kinematics) are tested. Two main hovering modes are investigated and compared with experimental and other computational efforts. The analysis shows that, while delayed-stall and rapid pitch-up mechanisms are responsible for most of the lift generation at a Reynolds numbers of O(100) and stroke amplitudes of O(1 chord), other mechanisms, including wake-capturing, are identified to contribute to the overall lift/drag force generation. The effect of the Reynolds number on hovering airfoil aerodynamics is also probed. xi

PAGE 12

CHAPTER 1 INTRODUCTION TO FLAPPING FLIGHT 1.1 Introduction From early times, the idea of utilizing the thrust generated from flapping wings for man-made flying vehicles emerged from observations of birds. One of the first attempts to explain and implement the mechanism of thrust generated from flapping wings to a flying vehicle was made by Leonardo da Vinci in 1490. Numerous attempts to develop flying vehicles based on flapping wing aerodynamics were made at the end of the 19th century and the beginning of the 20th century (Lilienthal, 1889; Lippisch, 1960), but most of them failed due to inadequate scientific knowledge. One of the first attempts to explain the physics of flapping wings was made by Knoller (1909) and Betz (1912). They noted that the wing oscillations produce both thrust and lift components of the aerodynamic force. In 1922, Katzmayr experimentally measured the thrust generated by a rigid wing in oscillating flow, confirming the Knoller-Betz effect. Based on Prandtls formulation (Prandtl, 1924), Birnbaum (1924) developed one of the first analytical approaches for calculating thrust and efficiencies of flat plates in plunging and pitching motion at moderate reduced frequencies. Karman and Burgers (1935) suggested an explanation for drag or thrust generation of oscillating airfoils based on the position and orientation of the vortices shed in the wake. During the second half of the 20 th century, research on flapping wing propulsion intensified, fueled by the fact that numerous biological organisms using flapping wing propulsion (in air or in water) possess superior aero-hydrodynamic properties and use 1

PAGE 13

2 energy more efficiently compared with those of the man-made flying vehicles, surface and underwater ships. Even with the rapid advance of aeronautical knowledge in the last 100 years, insects and birds still amaze us with their impressive performances. For example, a supersonic aircraft such as the SR-71 Blackbird flying nearly 2000 mph covers about 32 body lengths per second while a common pigeon flying at 50 mph covers nearly 75 body lengths (Shyy et al., 1999). Therefore, natures flapping wing systems are a precious source of inspiration for efficient propulsion and lift generation devices used in micro air vehicles. The use of flapping wing devices is very attractive, because they are able to operate in different regimes of motion, combine the function of control device, stabilizer and propulsor, and provide high maneuverability. However, the development of micro air vehicles by directly copying certain properties of flying organisms is unlikely to be successful. It is important to thoroughly analyze the applicability of the solutions found by nature for technical use, and this requires a multi-disciplinary effort involving aero-hydrodynamics, advanced materials, aero/hydro/servo/elasticity, micro-electro-mechanical systems (MEMS) technologies, etc. 1.2 Flapping Wing in Nature If we define energy cost as the energy required to move one kilogram a distance of one kilometer, and plot it against the body mass of birds, terrestrial animals, and fish (Schmidt-Nielsen, 1972; Goldsping, 1977), one can see that for a given body size, running is the most expensive in terms of energy costs, while flying and swimming are the cheapest. This may be the main reason why nearly a million insect species and more than 9,000 warm blooded animals have taken to the skies. To counter the gravity force

PAGE 14

3 and propel themselves against the aerodynamic drag, these creatures exclusively use the beating motion of wings. Different aspects of flapping wing flight for both birds and insects, such as type of flight, beating frequency, and scaling, are presented next. 1.2.1 Wing Motions Azuma (1992) identified four fundamental motions of a beating wing (Fig. 1-1): 1. an out-of-plane motion called flapping (flapping angle ); 2. an in-plane motion called lagging (lag or azimuthal angle ); 3. a twisting motion of the wing pitch called feathering (feathering angle ); 4. an alternatively extending and contracting motion of the wingspan called spanning. The stroke plane is defined by the elevation angle, s The feathering motion is also called supination for positive pitch and pronation for negative pitch. s Wing element Z X Y Horizontal axis Vertical axis s Figure 1-1. Wing beating motions. Sketched from Azuma (1992).

PAGE 15

4 Insects employ the same fundamental wing motions, except for spanning, since the structure of the insect wings do not allow it. 1.2.2 Flight Modes in Birds Birds can adopt different flight modes to fulfill their objectives in flight by adjusting the configuration of their wings or wing profile. This allows them to obtain the optimal wing configuration for each flight mode and each phase of wing beat (Azuma, 1992). In general birds use the following types of flight (Azuma, 1992): Cruising flight. Cruising flight is a steady level flight in which the almost fully extended wing acts to give lift force at the inner part of the wing and thrust at the outer part of the wing with minimum power consumption. Gliding or soaring flight Gliding or soaring flight is performed without beating the wings, except while maneuvering. The wings are fully extended in order to increase the wings aspect ratio and obtain a minimum gliding angle (max C L /C D ) or a minimum rate of descent (max ). Here C 3/2/LCC D L and C D are the lift and drag coefficient, respectively. Diving flight. Diving flight is a steep descent used for preying or landing. The wings are partially retracted to reduce drag. The birds control the speed and direction with the help of the tail. Hovering flight. Hovering flight is the flight mode in which birds can temporarily maintain their position in calm air. During hovering, the bird wings move in a complex pattern, with the wing tip describing an elliptical orbit. During the downstroke, the wings are extended to provide more lift, while on the upstroke the wings are flexed backward to reduce drag. This is usually called avian stroke (Azuma, 1992) and is illustrated in Fig. 1-2(A). Hummingbirds and insects can sustain longer periods of hovering flight by

PAGE 16

5 beating the wings almost horizontally with very high frequency, while the wing tip is moving in a figure 8 pattern. This is called insect stroke and is shown in Fig 1-2(B). Take-off flight This is the portion of the flight in which birds need to acquire enough speed to utilize the aerodynamic forces for a normal flight. Small birds can use their legs to jump in the initial stages of take-off, whereas large heavy birds take long runs against the wind to reach take-off speed. Landing. In this flight mode, birds need a large drag force to bleed their speed, while maintaining the lift necessary to sustain their weight at low speed. The wing shape during landing prevents flow separation at high angle of attack, and usually the tail is used as a flap to increase lift. 1.2.3 Beating Frequency and Scaling The mode and frequency of the beating motion depend on body size and shape, flight mode, and flight regime. Therefore, different species use very different beating motions. One thing they share is the ability to select the flapping mode and frequency that lead to optimal power consumption of the respective flight mode for normal flight. downstroke (power stroke) upstroke downstroke upstroke (recovery stroke) B) A) Figure 1-2. Hovering flight. A) Avian stroke. B) Figure pattern observed in hovering flight of hummingbirds and most of the insects. Sketched from Azuma (1992).

PAGE 17

6 The range of beating frequency and Reynolds number varies over a large range (as seen in Table 1-1), where the body mass, beating frequency and Reynolds number are presented for different species of birds and insects. In hovering flight, the maximum velocity of the wing tip is given by the formula (Azuma, 1992): 1 bfU (1-1) where f is the flapping frequency, b is the wing span, and 1 is the flapping amplitude. The maximum Reynolds number at the tip is defined as (Azuma, 1992): 121)/()/(/ReARfbcfbcU (1-2) where c is the mean chord defined by ARb/ c and AR is the wing aspect ratio. Then the reduced frequency can be defined as (Azuma, 1992): ARARUcfk/1)/(1/1 (1-3) The above equation states that a bird or an insect with a high aspect ratio wing flies at a small reduced frequency. Generally, larger flying creatures fly by gliding or slow beating, while smaller ones fly by strong beating at high frequency. Two main factors limit the beating frequency in flapping flight. One is the structural limitation of the bones (for birds), the supporting veins (for insects), or the muscles, wherein a higher load can lead to structural failures. Another constraint is the power available from the muscles. These limitations create the upper and lower bound for the wing beat frequency. The power consumed to overcome the inertial torque of the wing is proportional to the product of the moment of inertia and the cube of the beating frequency as shown below (Azuma, 1992): 35/3inPIfmf 3 (1-4)

PAGE 18

7 Table 1-1. Variety of beating frequency and Reynolds number for birds and insects. Species Mass m [kg] Beating frequency, f [Hz] Reynolds number, Re Chalcid wasp (Encarsia formosa) 2.5x10 -8 370 2.0x10 1 Fruit fly (Drosophila virilis) 2.0x10 -6 240 2.0x10 2 Fly (Musca domestica) 1.2x10 -5 190 6.4x10 2 Mosquito (Aedes nearcticus) 3.5x10 -6 320 1.3x10 3 Honeybee (Apis mellifica L.) 7.8x10 -5 250 2.3x10 3 Bumblebee (Bombus terrestris) 8.8x10 -4 156 4.0x10 3 Dragonfly (Anax parthenope) 7.9x10 -4 29 4.9x10 3 Locust (Schistocerca gregaria L.) 2.0x10 -3 20 1.0x10 4 Hummingbird (Patagona gigas) 2.0x10 -2 15 1.5x10 4 Sparrow 3.0x10 -2 13 1.0x10 5 Pigeon 3.5x10 -1 6 2.0x10 5 Stork 3.5 2 4.0x10 5 Data collected from Azuma (1992), Osbourne (1951) and Jensen (1956). Azuma (1992) states that the maximum power generated by a muscle during one cycle of flapping can be considered proportional to the product of the mass and the frequency: aPmf (1-5) Equating the maximum available power and the consumed power, the upper limit of the wingbeat frequency is given by: 5/331/3maxinaPPmfmffm (1-6) This same relation was specified by Hill (1950), who was one of the first to observe the direct proportionality between the mechanical power produced by a particular flight

PAGE 19

8 muscle and the contraction frequency. This has made flapping frequency an important parameter when trying to describe the theory behind flapping wings. To find the lower limit of wingbeat frequency, one can consider the force equilibrium in cruising flight as the flight mode with minimum power consumption, and consequently, with the lowest flapping frequency. In the vertical direction, the lift force cancels the weight and the following proportionalities can be written: 221/2LLWUSCmgUSmUmSUm 1/6 (1-7) where the geometric similarity for the wing area S was used () (Azuma, 1992). 2/3Sm In the horizontal direction, the thrust and drag are in equilibrium and the following proportionalities are obtained: 23211/()/2ADTDPUUSCSbfUUSfUbfUm 1/3 (1-8) where P A is the aerodynamic power and the geometric similarity for wingspan b was used (bm) (Azuma, 1992). Combining equations (1-7) and (1-8), the lower bound for wingbeat frequency is given as: 1/3 1/31/6fUmm (1-9) The same expression for lower bound of frequency was obtained by Shyy et al. (1999), starting with the assumption that the lower frequency is experienced by most of the birds in slow forward flight or hovering. Pennycuick (1996), using dimensional analysis, gave the following relation between the beating frequency f, wing span b, and the wing area S for birds in cruising flight: 3/823/241/33/81/2 f mbSg (1-10)

PAGE 20

9 where is the air density and g is the gravitational acceleration. Introducing geometric similarity expressions for the wing span and wing area (,) (Azuma, 1992), equation (1-10) leads to: which is consistent with the result given in equation (1-9). 1/3bm 2/3Sm 1/6fm Based on statistical data from birds and insects, Greenewalt (1962) showed that the wingbeat frequency, f, varies with the wing length l w (~ b/2) according to the relation: 1.15w f l = constant. For insects, the stroke frequency is widely scattered when plotted against the body mass but is confined between the and mline, as can be seen in Fig. 1-3. 2/1m 6/1 10-8 10-7 10-6 10-5 10-4 10-3 10-2 100 101 102 103 104 log(mass [kg])log(wingbeat frequency [Hz])Wing beat frequency vs body mass for different insect species different insect speciesf = m-1/2f = m-1/3f = m-1/6 Figure 1-3. Wing beat frequency versus the body mass for different insect species. (Data collected from Azuma, 1992). In his work, Azuma (1992) showed that the ratio between necessary power for flapping flight, P n and available power, P a increases with the mass of the flying animal,

PAGE 21

10 and therefore confirms the fact that the size of flying birds and insects is inevitably limited. A review of flapping wing kinematics and aerodynamic models for analyzing the lift, drag and power at low Reynolds numbers, as well as the role of flexible structures capable of adjusting to the free stream variations can be found in Shyy et al. (1999). 1.3 Insect Flight The main objective of this work is to develop a better understanding of the flapping wings aerodynamics at low Reynolds numbers, for possible applications to micro air vehicle design. Therefore, the main focus of this study is on insect flapping wings. First, the main differences between birds and insects flight regimes are presented. Second, the main efforts in experimental research of the insect flight mechanisms are briefly reviewed. The immense progress in computational methods and computer power allows us to perform simulations of the flow over an insect flapping wing, and get more insight into the unsteady flight mechanics of insect flight. These efforts are presented next, followed by a description of the main mechanisms responsible for the aerodynamic force generation. 1.3.1 General Characteristics Based on their sizes, flying insects operate over a range of Reynolds numbers from approximately 10 to 10 3 High beat frequencies and small sizes (0.25 mm to 0.1 m) result in low wing loadings (less than 0.1 N/m 2 ) and high reduced frequencies (more than 0.1). Even if the inertial forces are at least one order of magnitude larger than viscous forces, viscous effects have an important role in structuring the flow, and thus cannot be ignored. Birds and insects use the aerodynamic forces differently. Flying at high Reynolds numbers, birds only use the lift force, while reducing drag. Insects use both lift and drag,

PAGE 22

11 which are strongly coupled by unsteady effects generated by the low Reynolds number, high flapping frequency, and relatively low wing aspect ratio. Generally, insects show remarkable flight agility and maneuverability, hover for long periods, fly in any direction (even backwards), and are capable of landing upside down on leaves or flowers (Alexander, 1986; Azuma, 1988). For an airfoil, as the Reynolds number decreases, the maximum C L /C D ratio decreases and aerodynamic performances deteriorate. For birds, the feathers overlap and form a profile that resembles a smooth technical airfoil. Insect airfoils have rough surfaces such as well-defined cross-sectional corrugations (dragonfly) or scales on the wing surface (butterfly and moth). Nachtigall (1976) reported that the removal of the scales from a moth flying at a Reynolds number of 3x10 3 reduced lift by an average of 15% compared with the scaled wing. Kesel (2000) compared the aerodynamic characteristics of dragonfly wing sections with technical airfoils and flat plates. In her experimental study, she concluded that corrugated airfoils, such as those seen in dragonflies (Fig. 1-4), have very low drag coefficients closely resembling those of flat plates, while the lift coefficient was much higher than those of flat plates. The measurements were done at two Reynolds numbers: 7.8x10 3 and 1x10 4 The cross-sectional configurations vary greatly along the longitudinal axis of the wing, producing different aerodynamic characteristics. 1.3.2 Experimental Efforts in Studying the Insect Flight Mechanisms Because of their high wing beat frequencies and small size, it is quite difficult to measure the wing motions of free flying insects. Marey (1868) was the first to use slow-motion photography to track tethered insects wingtip motion. More recent methods have

PAGE 23

12 employed high-speed videography (Willmott and Ellington, 1997b), which offers greater light sensitivity. Measuring the aerodynamic forces history on insects is a very challenging task. At best, the forces have been measured on the body of a tethered insect, making it very difficult to separate the inertial forces from the aerodynamic forces (Cloupeau et al., 1979; Buckholz, 1981; Wilkin and Williams, 1993). Figure 1-4. Cross-sectional corrugations of dragonfly wings. Left wing is the hind wing and to the right is the fore wing. Picture taken by Dragos Viieru. Experiments on tethered insect flight give an important insight into the flow structures in the wake, but due to their dimensions the visualization results of the flow near the flapping wings are not clear enough, and in addition, tethering can alter the wing motion as compared with free flight conditions. Also, it is very difficult to measure the forces and moments generated by the insect wings in all degrees of freedom, and with enough temporal resolution required to link the wake structures with the aerodynamic force generation. To overcome these difficulties, many researchers have developed mechanical models of insect wings, matching the Reynolds number and reduced frequency parameter (body velocity/wing velocity) to that of an actual insect. These

PAGE 24

13 robotic devices have provided new insights into the vortex structures and aerodynamic forces of insect flight. A summary of the mechanical and robotic models is presented in Table 1-2. 1.3.3 Computational Efforts in Studying the Insect Flight Mechanisms The recent advances in computational methods and increasing computer power offer a new approach in solving the mechanism of insect flight. The computational methods applied to solving the fluid flow equations are more rigorous than the simplified analytical solutions, but still require empirical data for validation of wing kinematics. Nevertheless, based on more relevant input data from experiments, the computational fluid dynamics (CFD) models offer a better insight in the unsteady mechanisms of flapping flight. One of the first computational studies modeling the flight of a tethered hawkmoth was done by Smith (1996). In his investigation, he used an unsteady panel method for solving the governing equations and accounted for the wing flexibility using a finite element model. The aerodynamic forces at the center of the aerodynamic panels were transferred to the structural control points using a load transformation matrix computed using a finite surface spline method (Appa, 1989). The simulation results showed the effects of both the aerodynamic and inertial forces on flapping, flexible wings. The study also established the importance of including the wake in analysis of the unsteady aerodynamics of flexible wings undergoing arbitrary large-scale motion. Liu and Kawachi (1998) were the first to attempt a full unsteady Navier-Stokes simulation of the flow over a hawkmoth wing on a structured grid. To avoid the artificial mass fluxes due to grid movement, the Geometric Conservation Law (Thomas and

PAGE 25

14 Lombard, 1979) was implemented. The kinematic model of the flapping wing was based on the kinematic analysis of a hovering hawkmoth by Willmott and Ellington (1997a). Table 1-2. Experimental robotic and mechanic devices for insect flapping wing study. Author, Year Wing model based on: Main features Bennett (1970) Beetle (Melolontha vulgaris) 3-D mechanical model flow visualization only wing incidence can modified do not allow force measurements Spedding and Maxworthy (1986) Generic 2-D mechanical model for fling mechanism only wing rotation allowed allow measurements of the aerodynamic forces using a force transducer Saharon and Luttgers (1987), (1988) and Savage et al. (1979) Dragonfly 3-D mechanical model flow visualization Ellington et al. (1996) Hawkmoth 3-D robotic model good flow visualization (model wing =~10 times insect wing) computer-controlled wing kinematic do not allow force measurements Dickinson et al. (1999) and Lehmann (2000) Fruit fly 3-D robotic model good flow visualization (scaled-up wing) computer-controlled wing kinematic allow measurements of shear forces using a 2-D force transducer To facilitate the resolution of viscous flow at the leading edge, trailing edge and wing tip, an O-O type grid topology was used. This study confirmed the smoke streak patterns observed in both real and dynamically scaled model insect flight (Ellington et al., 1996), and accurately predicted the complex vortex structure and the importance of the span-wise flow in stabilizing the spiral leading-edge vortex as a lift-enhancement mechanism, as one can see in Fig. 1-5. Ramamurti and Sandberg (2002) studied the unsteady flow over a rigid fruit fly wing based on a dynamically scaled model using a finite element solver and a re-meshing technique. A detailed description of the mesh movement algorithms is given in Ramamurti et al. (1994). The effect of phasing between the translational and rotational motions was studied.

PAGE 26

15 To p View Side View Tip vortex LEV Side View A B Figure 1-5. Leading edge vortex on a hawkmoth wing. A) Numerical result. B) Smoke-visualized flow during downstroke. Reprinted from Liu et al. (1998) with permission from The Company of Biologists Ltd. It was observed that the highest peak in the thrust force is obtained when the wing rotation is advanced with respect to the stroke reversal, followed by the case in which the wing rotation is symmetrical with respect to the stroke reversal. The lowest thrust is obtained when the wing rotation is delayed. The computed thrust and drag roughly matched the experimental results measured on a dynamically scaled model (Dickinson et al., 1999), and showed the importance of the rotational mechanisms in accurately describing the force time histories. Another analysis on the unsteady aerodynamics of a model fruit fly wing, including the energy consumption aspect, was conducted by Sun and Tang (2002a, 2002b) using an implicit algorithm that is second-order accurate in time and space on a structured grid.

PAGE 27

16 The wing movement was handled using a time dependent coordinate transformation between an inertial reference frame and the non-inertial, body fixed frame of reference. Three mechanisms were observed to account for the large lift: the rapid acceleration of the wing at the beginning of a stroke, the delayed stall during the stroke, and the fast pitching-up rotation of the wing near the end of the stroke (advanced rotation). The computed unsteady lift coefficients were smaller than those measured on a robotic wing by Dickinson et al. (1999). This was partly explained by using a wing with a smaller aspect ratio than that of the robotic wing. However, the overall computed lift coefficient pattern was very similar to the experimental results. Recently, Gilmanov and Sotiropoulos (2005) performed simulations of flow over the fruit fly (Drosophila). The wing kinematics and flow parameters were those considered in the experimental work of Birch and Dickinson (2003). Their results showed a reasonable agreement between the computed force history and the experimental values (Birch and Dickinson, 2003). They numerically solved the three-dimensional, unsteady, incompressible Navier-Stokes equations in a Cartesian domain using the immersed boundary method, which is capable of handling complex geometries moving with prescribed kinematics. The surface of the body moving in a fixed Cartesian domain is discretized using an unstructured triangular mesh. The nodes of the surface mesh constitute a set of Lagrangian control points used to track the movement of the body. At each time step the influence of the body on the flow is accounted for by applying boundary conditions at Cartesian grid nodes located in the exterior, but in the immediate vicinity of the body by reconstructing the solution along the local normal to the body surface. To discretize the governing equations, second-order finite difference formulas

PAGE 28

17 are used on a hybrid staggered/non-staggered grid layout. The discrete equations are integrated in time using a second-order, dual-time-stepping, artificial compressibility iteration scheme. The aerodynamic force generation and mechanical power requirements of a dragonfly in hovering flight were studied by Sun and Lan (2004). Two O-H type grids were generated around the forewing and the hindwing, and they moved in a fixed Cartesian grid. The wing kinematics were based on experimental data of a hovering dragonfly (Norberg, 1975; Reavis and Luttges, 1988). The computational results showed two large vertical force peaks in one flapping cycle generated by a new vortex ring containing downward momentum forming on the forewing or the hindwing during each downstroke. The interaction between the forewing and hindwing was not found to be very strong, and even detrimental to the vertical force generation. The study also showed that the dragonfly uses drag as a major source of vertical force (approximately 65% of the total vertical force is contributed by the drag) during hovering. Similar results were obtained by Mittal et al. (2002) in their two-dimensional simulation of a hovering dragonfly. The solver employed a second-order accurate central difference scheme for spatial discretization and a mixed explicit-implicit fractional step scheme for time marching (Udayakumar et al., 1999; Ye et al., 1999; Udayakumar et al., 2001). Different values of phase lag between the forewing and hindwing were used. The results showed that, for the set of simulations studied, the overall efficiency of the paired wing system was lower than a single wing for all values of phase lag. For the single wing, the simulation indicated that the highest thrust efficiency was accompanied by the formation of an inverse Karman vortex street.

PAGE 29

18 Despite the importance of 3-D effects, comparison of experiments and computations in 2-D have also provided important insight. Wang (2000a, 2000b) demonstrated that the force dynamics of 2-D wings might be sufficient to explain some of the enhanced lift coefficients measured in insects. A summary of the main numerical studies of the insect flapping flight and governing equations solved are presented in Table 1-3. Table 1-3. Principal numerical studies of insect flapping flight Author, Year Fluid governing equations Method Moving boundary tracking Smith, (1996) 3-D potential flow Panel method Time depended boundary conditions Liu and Kawachi, (1998) 3-D Navier-Stokes in body-fitted coordinates Artificial compressibility (Chorin, 1968) Grid re-meshing. Analytical based on the initial grid and wing kinematics Wang, (2000a, 2000b) 2-D Navier-Stokes for vorticity in elliptic coordinates. Finite difference Time depending boundary conditions Mittal et al., (2002) 2-D Navier-Stokes equations Finite volume Immersed Boundary Method on fixed Cartesian grid Ramamurti and Sandberg, (2002) 3-D Navier-Stokes in ALE formulation (Arbitrary Lagrangian Eulerian) Finite Element Grid re-meshing Spring analogy and smoothing on unstructured grid Sun and Tang, (2002a, 2002b) 3-D Navier-Stokes in body-fitted coordinates. Artificial compressibility (Rogers et al., 1991; Rogers and Kwak, 1990) Time dependent coordinate transformation Sun and Lan, (2004) 3-D Navier-Stokes in body-fitted coordinates. Artificial compressibility (Rogers et al., 1991; Rogers and Kwak, 1990) Overset moving grid Gilmanov and Sotiropoulos, (2005) 3-D Navier-Stokes Finite difference Immersed Boundary Method on fixed Cartesian grid Mittal et al., (2006) 3-D Navier-Stokes Finite difference Immersed Boundary Method on fixed Cartesian grid 1.4 Unsteady mechanisms identified to enhance lift in insects 1.4.1 Analytical and quasi-steady models of insect flight Early models of insect flight analyzed the far-field wakes and could not be used to calculate instantaneous forces on airfoils, but could characterize the average forces. The

PAGE 30

19 most notable analytical models are the vortex models of Ellington (1978, 1980, 1984e) and Rayner (1979a, 1979b) that are applicable for hovering and forward flight. These vortex theories provide a method to compute the mean lift and induced power for a given circulation profile during the wingbeat. They consider the flapping wings as an actuator disc (pulsed actuator disc) that pushes air downwards at a certain rate. The rate of change of momentum flux within the downward jet must be equal to the insect weight, and therefore the circulation in the wake required to maintain the force balance can be calculated, leading to an estimation of the mean lift force. However, the momentum jet theory touches the insect flight on a superficial level, ignoring variations in kinematics and lift generation. A detailed description of these theories is presented in Rayner (1979a, 1979b) and Ellington (1984e). Azuma (1985) and Azuma and Watanabe (1988) used the local circulation method to study dragonfly and damselfly flight. This method incorporates the temporal and spatial (span-wise) changes in induced velocity and corrects the circulation due to wake. A more recent analytical model for flapping flight, assuming rigid wings, is presented in Zbikowski (2002). To predict aerodynamic forces from the motion of the wing, scientists developed simplified models based on the quasi-steady approximation. The main assumption is that the instantaneous aerodynamic forces on the flapping wing are equal to those generated during steady motion of the wing at identical instantaneous velocities and angles of attack (Ellington, 1984a). Therefore, the time dependence of the aerodynamic forces arises from the kinematics time dependence, and not from changes in the fluid flow. Also, the lift and drag coefficients are time invariant and derived typically from experiments in which

PAGE 31

20 model wings are tested over the range of angles of attack and flow velocities encountered during flapping motion. The theory was used with some success by Osborne (1951) and Jensen (1956) in their study of the desert locusts forward flight. However, the existing quasi-steady theory predicts lower lift coefficients than the required average lift coefficients for hovering as argued by Ellington (1984a, 1984b, 1984c, 1984d, 1984e) in a comprehensive review of the insect flight literature and based on a wide range of data available at that time. Direct force measurements in flying insects also validated his findings, indicating that insects do not use the steady-state principle (Ellington, 1984a, 1984b, 1984c, 1984d, 1984e). Based on experimental, computational and analytical models, some important unsteady mechanisms that play an important role in generating aerodynamic forces in insect flight can be identified, and they are described below. 1.4.2 Wagner Effect The Wagner effect is manifested when an inclined wing starts impulsively from rest. In this case the circulation around the wing does not immediately attain its steady state value. This mechanism that delays the development of the maximum circulation was proposed by Wagner (1925) and studied experimentally by Walker (1931). Based on Walkers results, the steady circulation is established after approximately 4 to 5 chord lengths. Ellington (1984d) defined the ratio of the mean distance traveled by the wing during a half-stroke to the mean chord as a measure of the number of chord lengths traveled. For the most insects this value is between 1.9 and 3.8; this fact indicates that steady circulation is never reached (Ellington, 1984a). It seems that, unlike other unsteady mechanisms, the Wagner effect would attenuate the unsteady forces below levels predicted by quasi-steady models. Studies done on 2-D wings by Dickinson and Gtz (1993) indicate that the Wagner effect is not very strong in the range of the

PAGE 32

21 Reynolds numbers experienced in insect flight, and most recent models of insect flight neglected the Wagner effect. 1.4.3 Leading-edge Vortex and Delayed Stall For thin airfoils traveling at high angle of attack the local viscous forces within the fluid near the leading edge are smaller than the pressure forces generated by the high fluid velocity, and the flow separates from the upper surface of the wing. The low pressure region behind the leading edge forces the flow to curl back and reattach on the upper wing surface behind the leading edge, leading to the formation of a leading edge vortex. In this case a suction force develops perpendicular to the wing chord, consequently enhancing the lift as well as the drag, as shown in Fig. 1-6. For a 2-D motion of the wing at high angle of attack, the leading edge vortex grows in size until reattachment is no longer possible. The leading edge vortex sheds into the wake, and as a result the lift drops and the wing stalls. But for several chord distances before stall, the attached leading edge vortex generates very high lift coefficients (delayed stall). At low Reynolds numbers, as the leading edge vortex is shed, another vortex is formed at the trailing edge that grows until it is also shed into the wake and a new leading edge vortex is forming. This dynamic process creates a wake of regularly spaced counter-rotating vortices known as the Karman vortex street (Dickinson and Gtz, 1993). One of the first visualizations of the leading edge vortex and the delayed stall mechanism in insect flight was done by Maxworthy (1979) on the model of a flinging wing. Since most insects flap their wings at high angle of attack, the leading edge vortex has an important role in lift generation.

PAGE 33

22 Lift Drag FSuction FResultant = FSuction + FNormal Figure 1-6. Flow around a thin airfoil. The presence of the leading edge vortex generates a suction force normal to the surface of the airfoil due to the low pressure zone associated with the vortex core. This force enhances the resultant force normal to the airfoil. The filled circle indicates the leading edge. Sketched from Sane (2003). In contrast with 2-D models, the leading edge vortex in 3-D flows was not shed even after many chord lengths and no Karman vortex street is present, as observed by Ellington et al. (1996) in his smoke visualization of the flow around both real and robotic hawkmoth wings at a Reynolds number in the range of 10 3 The mechanisms of 3-D leading edge vortex stability are still unclear. Ellington and co-workers observed a steady span-wise flow from the wings hinge to approximately three-quarters of the span towards the tip where the conical leading edge vortex detaches from the wing surface. This phenomenon was also observed by Liu and Kawachi (1998) in their computational study of the hawkmoth wing. The axial flow role in stabilizing the leading edge vortex was also observed in the case of back-sweep delta wings. However, the importance of axial flow in the stability of the leading edge vortex was questioned by experiments using a 3-D model wing based on the fruit fly flapping at a Reynolds number in the order of 100 (Birch and Dickinson, 2001). The measurements showed that the span-wise flow

PAGE 34

23 within the leading edge vortex core was quite small (2-5% of the average tip velocity), while a large axial flow on the wings upper surface behind the leading edge vortex was observed, that rolls into a large tip vortex. Recently, in an extended study of the effect of Reynolds number on the leading-edge vortex at mid-downstroke in hovering, Liu et al. (2005) found that the leading-edge vortices show significant variations in their structures as the Reynolds number changes. Figure 1-7 shows the streamline structures at three Reynolds numbers, corresponding to a hawkmoth in Fig. 1-7(A), a fruit fly in Fig. 1-7(B), and a thrips in Fig. 1-7(C). At a Reynolds number of 6,000 (corresponding to a moth) as observed in the previous studies (Ellington, 1996; Liu and Kawachi, 1998), an intense, conical leading-edge vortex is observed on the paired wings with a steady span-wise flow at the vortex core, breaking down at approximately three-quarters of the span towards the tip. At a Reynolds number of 120, corresponding to a fly (Fig. 1-7(B)), the break-down apparently disappears, and a leading-edge vortex is found to be connected to the tip vortex, forming a longer, leading-edge tip vortex. The vortex shows strong stability, while the span-wise flow at the vortex core turns to become much weaker compared with those at high Reynolds numbers, which is in agreement with the measurements of Birch and Dickinson (2001). Further reducing the Reynolds number down to 10 for the thrips case, a vortex ring connecting the leading-edge vortex, the tip vortex and the trailing vortex is observed (Fig. 1-7(C)). Their 3D structure looks more like a cylinder than the conical form seen at high Reynolds numbers. Inspecting the momentum equation, one can see that the pressure-gradient, the centrifugal force, and the coriolis forces are all together responsible for the leading-edge

PAGE 35

24 vortex stability. Their roles at different Reynolds numbers should be studied to help shed light on these findings. A) B) C) Figure 1-7. Leading edge vortex (LEV) structures at different Reynolds numbers. A) Reynolds number = 6000. B) Reynolds number = 120. C) Reynolds number = 10. Adopted from Liu et al. (2005). In another effort, Sane and Dickinson (2002), by placing a fence on the upper surface of a robotic model wing, showed that the leading edge vortex (LEV) still exists on the wing. Based on this observation, they suggested that the spanwise flow should not be responsible for the LEV on the upper wing surface. 1.4.4 Rotational Forces When a wing rotates about a span-wise axis and translates at the same time, it generates a rotational circulation in the fluid to counteract the effects of the rotation. The additional circulation is proportional to the angular velocity of the rotation and generates rotational forces, which, depending on the direction of the rotation, can be added or subtracted from the net force due to translation. Insects rotate their wings around the span-wise axis near the end of every stroke in order to obtain a positive angle of attack for the next half-stroke and thus, the forces generated during rotation are very important. Direct force measurements on tethered fruit flies have shown that maximum forces are produced during wing rotation, suggesting that the rotational forces generated during stroke reversal represent a significant part of the total lift (Dickinson and Gtz, 1996).

PAGE 36

25 Similar rotational force peaks were observed in CFD simulations by Sun and Tang (2002a). Both CFD simulations and experimental force measurements indicate that the forces during flight are sensitive to the phase of wing rotation. Advancement in wing rotation at the end of a half-stroke generates the maximum mean lift coefficient, while delaying the rotation reduces the mean lift coefficient by almost 42% (Dickinson et al., 1999). Dickinson et al. (1993) suggested that rotation timing is used by insects during steering maneuvers; the insect advances the rotation on the wing outside the turn, and delays it on the wing inside the turn. 1.4.5 Wing-wake Interaction In contrast with wings of aircraft that move through still air, the wings of hovering insects intercept the wake created by the wings own movement. The interaction of insect wings with the shed vortex of prior strokes was observed with both force measurements and flow visualization on a robotic model based on a fruit fly (Dickinson et al., 1999). During stroke reversal both leading and trailing edge vortices are shed, and consequently a strong velocity field is generated between these two vortices. As the wing reverses direction, it encounters the enhanced velocity field that transfers momentum to the wing, thereby enhancing the maximum force during stroke reversal. This mechanism is illustrated in Fig. 1-8. The magnitude and relative strengths of shed vortices are strongly dependent on the wing kinematics immediately before and after stroke reversal. The model wing produces positive wake capturing lift when the wing rotation precedes stroke reversal. Similar force peaks were observed by Sun and Tang (2002a) in their CFD simulations of a fruit fly based on kinematics similar to those in Dickinson et al. (1999), but they attribute these peaks to the fast pitching-up rotation of the wing, and not to the wing-wake interaction. They also provided some experimental cases to prove

PAGE 37

26 their conclusion. Later, Walker (2002), using a different CFD method, reached the same conclusion. U A B C D Transferred momentu m Local flow Downwas E F Figure 1-8. Momentum transfer due to wing-wake interaction. A) Wing is in steady translation phase. B) Leading and trailing edge vortices are generated as the wing rotates around a spanwise axis. C) Velocity field induced by the shed vortices. D) Wing rotates and starts the reverse stroke. E) Wing reverses the stroke and encounters the induced velocity field. Fluid momentum is transferred to the wing that generates a peak in the aerodynamic. F) Wing is in reversed stroke steady translation. The filled circle indicates the leading edge. Sketched from Sane (2003). The wake capturing mechanism may explain the insects ability to extract energy from its own wake, recovering energy from air lost during the previous stroke and therefore improving the overall force production efficiency. Srygley et al. (2002) observed, during a smoke visualization experiment, the extraction of kinetic energy from the wake behind freely flying butterflies during take-off. However, the contribution of the wing-wake interaction mechanism to the peak force generation during wing rotation is still not clear, and requires further investigation.

PAGE 38

27 1.4.6 Clap and Fling Mechanism The clap-and-fling mechanism is the physical interaction of the left and right wing during dorsal stroke reversal. This mechanism was proposed by Weis-Fough (1973) to explain the lift generation in the chalcid wasp. Its schematic is presented in Fig. 1-9. Figure 1-9. Schematic of clap and fling mechanism. Based on Weis-Fough (1973). This mechanism was also found in fruit flies (Gtz, 1987), butterflies (Brodsky, 1991), and locusts (Cooter and Baker, 1977). During the fling phase, preceding the down-stroke, the fluid is pulled into the gap formed by the moving wings, generating strong leading edge vortices while the development of the trailing edge vortices is stopped by the contact of the wings trailing edges. When the wings translate away from each other, the low pressure region between them accelerates the surrounding fluid, providing a build-up of circulation. A detailed theoretical analysis of the clap-and-fling mechanism can be found in Lighthill (1973) and Sunada et al. (1993), and experimental observations can be found in Bennett (1977), Maxworthy (1979), and Spedding & Maxworthy (1986). The overall benefit of the clap and fling remains uncertain in many studies of insect wing kinematics. Recently, the unsteady mechanisms of lift enhancement in insect flight were reviewed in Lehmann (2004). A review of the current models as well as of the unsteady

PAGE 39

28 mechanisms of insect flight is presented in Sane (2003). Recent findings in the unsteady aerodynamics of insect flight at intermediate Reynolds numbers (10-10 4 ), with focus on two-dimensional unsteady mechanisms in uniform and accelerated motions, forward and hovering flight can be found in Wang (2005). A comprehensive survey of the research of flapping wing propulsors and their application to man-made vehicles has been presented by Rozhdestvensky and Ryzhov (2003). 1.5 Objective and Outline The main objective of this work is to investigate the aerodynamic force generation mechanisms for fixed and flapping wings micro air vehicles at low Reynolds numbers. In the light of the insect flight analysis presented above, it is clear that a better understanding of the unsteady mechanisms for lift generation is very important if we want to explore the development of micro and nano air vehicles. Therefore, in this work, the focus is to better understand biological lift generation by flapping wings at low Reynolds number, as well as to point out some of the characteristics of small micro air vehicles with fixed rigid or flexible wings. In this chapter, the flapping wing in nature was discussed with a focus on insect flight. Main experimental and computational research efforts were presented, followed by a presentation of the main unsteady mechanisms identified so far that are responsible for lift generation in insects. In chapter 2, the numerical framework for solving fluid-structure interaction problems with moving boundaries is presented. The fluids governing equations are presented as well as the numerical implementation of those equations in a pressure-based fluid solver. The finite element structural solver and the governing equations for solids are also presented. The fluid-structure interaction problem leads to deformation in the

PAGE 40

29 computational grid and therefore algorithms to automatically take into account the grid movement are described. Since the fluid and structure solvers require different grid types, an interpolation technique between the structured CFD grid and the unstructured FEM (finite element method) mesh is mentioned. Finally, a generic algorithm for solving the fluid-structure interaction problem is presented. In chapter 3, the flow over fixed rigid and flexible wing MAVs is analyzed. The low aspect-ratio wings used in the University of Florida MAV design generates strong tip vortices that have an important impact on the wings aerodynamic performances. The wing tip vortex effect is studied along with means to reduce the induced drag and increase the rigid wing efficiency. Lastly, the aerodynamic performances of a flexible wing are compared with those of a similar rigid wing. The two-dimensional hovering of an elliptic airfoil is investigated in chapter 4. The results are compared with other computational work to assess the capabilities of the fluid solver to handle this type of problem. Also, basic unsteady mechanisms for lift generation are explored, such as delayed-stall and wake-capturing mechanisms. Finally we conclude with the summary and the directions for the future work in chapter 5.

PAGE 41

CHAPTER 2 COMPUTATIONAL METHODS FOR COUPLED FLUID-STRUCTURE INTERACTIONS 2.1 Overview This study investigates the aerodynamics of fixed-wing MAVs with both rigid and flexible wings, as well as the insect flapping flight under low Reynolds number conditions. During flight, the flexible wing changes its shape as a response to the external forces and, at the same time, its shape variation affects the pattern of its surrounding fluid flow. Under steady freestream the membrane wing exhibits high frequency vibration, observed both experimentally (Waszak et al., 2001) and numerically (Lian et al., 2003c; Lian and Shyy, 2005). Also, the insect flapping flight is dominated by unsteady mechanisms such as vortex shedding, dynamic stall, and wing-wake interaction. It is desirable to have a computational capability to accurately capture the transient complex behavior of fluid flow around flexible and flapping wings. Numerous approaches have been proposed for solving unsteady flows. Pulliam (1993) incorporated subiterations to improve the time accuracy of conventional implicit schemes. Jameson (1991) and Rumsey et al. (1996) developed a dual time stepping method used in the multigrid framework to improve the efficiency and accuracy of traditional factored schemes. Issa (1985), and Oliveira and Issa (2001) proposed the Pressure-Implicit Splitting of Operators (PISO) method to improve convergence. For a review of these methods and some other frequently adopted approaches, see Shyy and Mittal (1998). 30

PAGE 42

31 In this work, the three-dimensional Navier-Stokes equations for incompressible flow written in curvilinear coordinates are solved. Two popular methods to solve the Navier-Stokes equations for incompressible flow are employed. One is the SIMPLEC method (van Doormaal and Raithby, 1985) that is an enhancement of the SIMPLE algorithm (Patankar and Spalding, 1972; Shyy, 1994), and the other one is the PISO method (Issa, 1985). Both methods belong to the pressure-based family (Shyy, 1994) that uses an artificially derived pressure or pressure correction equation. In this chapter, the fluid solvers adopted for moving boundary problems are briefly introduced, as well as the importance of the boundary conditions and the geometric conservation law (Thomas and Lombard, 1979). Different moving grid techniques used to automatically re-mesh the grid are investigated. Finally, the membrane structural solver and the coupling between the fluid solver and structural solver are discussed. 2.2 Governing Equations for Fluid Flows The three-dimensional Navier-Stokes equations for incompressible flow in Cartesian coordinates can be written, using indicial notation, as follows: 0jjuxt (2-1) jjiiijjixxpuuxut (2-2) where x i is the position vector, t is time, is density, u i is the velocity vector, p is pressure, and ij is the viscous stress tensor. The constitutive relation between stress and strain rate for a Newtonian fluid is used to link the components of the stress tensor to velocity gradients:

PAGE 43

32 ijllijjiijxuxuxu32 (2-3) where is the molecular viscosity. For arbitrary shaped geometries, the Navier-Stokes equations are transformed into generalized curvilinear coordinates (, ), where ),,(,,,zyxzyx and ),,(zyx The transformation of the physical domain (x, y, z) to the computational domain (, ) is achieved by the following relations (Tannehill et al., 1997): 3332312322211312111fffffffffJzyxzyxzyx (2-4) where f ij are the metrics terms and J is the determinant of the Jacobian transformation matrix given by: ,,detdet,, x xxxyzJyzzz yy (2-5) To solve for the Navier-Stokes equations in curvilinear coordinates, the finite-volume formulation is adopted. In this approach, both Cartesian velocity (as primary variables) and contravariant velocity components are employed. The contravariant velocities are used to evaluate the flux at the cell faces and to enforce the mass continuity in the pressure-correction equation. The expressions for the metrics, the determinant of the Jacobian transformation matrix, and the fluxes at the cell faces for the three-dimensional Navier-Stokes equations in the generalized body-fitted curvilinear coordinate system (, ) are given in Thakur et al. (2002) and Shyy (1994). When the

PAGE 44

33 governing equations are considered under a moving grid framework, the grid velocities should be included in the flux computations as described in Shyy et al. (1996). 2.3 Fluid Flow Solvers The SIMPLEC method (van Doormaal and Raithby, 1985) used in this study to solve the governing equations is an enhanced variant of the SIMPLE algorithm initially proposed in Cartesian coordinates by Patankar and Spalding (1972). The extension to body-fitted curvilinear coordinates capable of handling complex geometry computation was discussed by Shyy (1994). An attractive feature of these methods is that the solution procedure is sequential in nature and can therefore easily accommodate a varying number of equations depending on the physics of the problem involved. The implementation of the method employs a control volume approach and uses a non-staggered arrangement for the velocity components and the scalar variables (i.e. pressure). The Cartesian velocity components are computed from the respective momentum equations. The cell fluxes and pressure fields are corrected using a pressure correction equation, which is derived by manipulating the continuity and momentum equations. The iterative correction procedure leads to a divergence-free velocity field within a desired convergence tolerance, therefore enforcing the pressure-velocity coupling. The SIMPLEC method uses coordinated under-relaxations for the momentum and pressure corrections to improve the convergence speed (intrinsically slow in the original SIMPLE method). The other method used in this work is the PISO method. Contrary to the iterative SIMPLE family of methods, the PISO method uses a non-iterative process. Here the pressure-velocity coupling is handled by splitting the solution process into a series of predictor and corrector steps. The PISO method is generally more efficient for transient flow computations (Issa, 1985; Thakur et al., 2002).

PAGE 45

34 To solve the Navier-Stokes equations, proper boundary conditions are required. At the interface between the fluid and solid structure, the no-slip condition is applied, which requires that the fluid velocity at the wing surface matches the surface velocity of the structure. These velocity boundary conditions are provided by computing the displacements of the solid structure. 2.4 Geometric Conservation Law To solve the governing equations in body-fitted curvilinear coordinates, a transformation matrix is used to facilitate the mapping of a physical flow region (x, y, z) onto a computational domain (, ). The Jacobian transformation matrix is defined as in equation (2-5). The determinant of the Jacobian transformation matrix represents the volume of the element in the transformed coordinate. In moving grid problems, the computational grid is changing with time and consequently the determinant of the Jacobian matrix, J, needs to be updated. Special procedures are required to compute the effective value of J at each time step; otherwise errors arise due to an inconsistent numerical implementation that would lead to the generation of artificial mass sources. As suggested by Thomas and Lombard (1979), in the process of updating the Jacobian (J) the following conservation law (GCL) needs to be satisfied: sVVdJdddWddddt (2-6) where V is the volume bounded by the closed surface S, and W s is the local velocity of the moving boundary surface S. Thomas and Lombard (1979) proposed an expression to evaluate the Jacobian (J) from the continuity equation for a constant density, uniform velocity field under a time dependent coordinate transformation, while maintaining the geometric conservation law:

PAGE 46

35 tJtJtJtJ (2-7) Integrating equation (2-7) with a first-order, implicit time integration scheme over the same control volume used for mass conservation leads to a finite volume discrete form of the above equation. This is used to update the Jacobian in a manner that respects the basic requirement of the geometric conservation in the discrete form of the conservation law when the grid is time dependent. Implications of the geometric conservation law on moving boundary problems were discuss by Shyy et al. (1996). 2.5 Moving Grid Technique 2.5.1 Master-slave Concept For moving boundary problems where a solid boundary (i.e. wing) moves inside a computational domain based on known kinematics (i.e. rigid flapping wing) or as a response of the structure to the flow around it (i.e. membrane wing), the grid needs to be adjusted dynamically during computation. To facilitate this, a moving grid technique needs to be employed. The actual process of generating a grid is a complicated task, so an automatic and fast algorithm to frequently upgrade the grid is essential. It is desirable to have an automatic re-meshing algorithm to ensure that the dynamically moving grid retains the quality of the initial grid, and avoids problems such as crossover of the grid lines, crossed cell faces, or negative volumes at block interfaces (in case of the multi-block grids). The moving grid technique can also be used for shape optimization studies (Lian et al., 2003b). Several approaches have been suggested to treat grid redistributions for moving grid computations. Schuster et al. (1990) used an algebraic shearing method in their study. The displacement of the moving surface is redistributed along the grid line which

PAGE 47

36 connects the moving surface to the outer boundary. This simple method gives good results for modest displacements and single block grids. For multiblock grid arrangements extensive user intervention is required. A robust method that can handle large deformations is the spring analogy method, which was first introduced by Batina (1989) for unstructured grids and later extended to structured grids by Robinson et al. (1991). In this method all edges of a cell, as well as the diagonals, are replaced by linear springs, each with the stiffness inversely proportional to a power p of the length of the connecting edges. Using the power p, one can control the stiffness of the spring and consequently control the amount of movement to avoid excessive mesh distortions. The iterative process necessary to find the displacement of all the internal points increases the computational time of this method, especially for large grids. The direct transfinite interpolation method was introduced by Eriksson (1982), and can generate grids for complex geometry. The method defines an interpolation function given known values on constant planes and function derivatives in the out of surface direction on the boundaries. The method is fast and efficient for structured grids, but the quality of the initial grid may not be always preserved, especially far from the boundaries. Hartwich and Agrawal (1997) propose the Master/Slave concept to expedite the grid regeneration process and preserve the grid continuity at the multiblock grid interfaces. In this study a moving grid technique proposed by Lian (2003a) and Lian et al. (2003d) is used to re-mesh the multiblock structured grid for fluid-structure interaction problems. For multiblock structured grids, for simplicity, CFD solvers often require point-matched grid block interfaces. This method is based on the master-slave concept

PAGE 48

37 and maintains a point-matched grid block interface while maintaining grid quality and preventing potential grid cross-over. When an object changes its shape, the master points, which are located at the moving body surface, move first, and then affect the distribution of the off-body points. One difficulty for a multiblock grid resides in the way in which the vertices of each block are moved, if a point-to-point match between two abutting blocks is required without overlap. For identical interfaces between two abutting blocks, such as the interface between block 2 and block 3 in Fig. 2-1(A), the edge and interior points can be obtained by a 3-stage interpolation once the corner vertices are determined. However, when the abutting blocks do not have an identical interface, such as of block 1 and block 2 in Fig. 2-1(A), the interpolation can cause a discontinuity at the interface. To avoid creating undesirable grid discontinuities, the off-body vertices of a grid block are linked to a surface point and thus move in a similar way. Therefore, for each off-body vertex (slave point), the nearest body surface point is defined as its master point. The distance between the slave point and its master is given by: 222msmsmszzyyxxr (2-8) where the subscript s represents a slave point, and m a master point. Once a slave point has its master point identified, the slave point moves according to the influence from its master. A simple formula to define the displacement is shown here: mmssxxxx ~ ~ (2-9) The tilde (~) indicates the new position, and is a decay function used to control and scale the movement to avoid a lower surface crossing over its opposite face when the displacement is large.

PAGE 49

38 Figure 2-1. Illustration of the moving grid technique. A) Initial geometry. B) Side view of the initial geometry. C) Modified geometry. D) Grid points distribution of the new shape. E) Same test but with higher stiffness. Adopted from Lian et al. (2003c). The decay function is a Gaussian distribution function given by: dmdv /,500minexp (2-10) where is a small number necessary to avoid division by zero when dm = 0, and .~~~,222222mmmmmmmvmvmvzzyyxxdmzzyyxxdv (2-11)

PAGE 50

39 The coefficient controls the stiffness of the movement. A larger value of causes the block to behave more like a rigid body as one can see in Fig. 2-1(E). More details about master and slave concept and its applications can be found in Hartwich and Agrawal (1997), Lian (2003a), and Lian et al. (2003c). Figure 2-2. Master-slave algorithm performance for large rotations The master-slave algorithm is highly automated, capable of handling large translational deformations while maintaining overall grid quality and allowing large translational displacements. The propagation distance of the moving wall perturbation is controlled by the spring stiffness coefficient. However, in the current formulation, large rotational deformations cannot be handled properly since no information is provided about cell skewness as the perturbation is propagated along grid lines. As one can see in Fig. 2-2, for large deformations, the grid quality deteriorates, leading to negative cell volumes. 2.5.2 Grid Re-arrangement Algorithm with Weighting Functions The flapping wing problem requires both large rotations and displacements of the airfoil, and an algorithm that can handle this type of problem needs to be developed. The algorithm proposed dissipates the shear stress created by the movement of the body

PAGE 51

40 toward the outer boundaries. A weight function is computed for each grid node based on a weighted average of the neighboring node positions. The main steps of the grid-rearrangement algorithm are as follows: i) Assign the boundary conditions. Assign weight function () = 1 (rigid rotation/translation) on the body surface and 0 on the outer boundary (no displacement) ii) Start the iterative process (k=1) iii) Compute the weight function value: 1,1/()1/()NnbrpmmmijNnbrpmmll where Nnbr is the number of neighboring nodes (8 for 2D see Fig. 2-3), l is the distance between the (i,j) node and the neighboring nodes, and p is a stiffness parameter (default = 2, but can vary). iv) check convergence criteria: 1,,kkijij < tolerance. If YES => new weight function, if NO => k=k+1 v) increment the new node position: xx 1,,newoldijijijdx (translation) => 1,,newnewijijijxxdyd (rotation), where dx, dy = translation increment, d = rotation increment. 8 7 6 5 4 3 2 m=1 i,j Figure 2-3. Current node neighbors position for grid re-arrangement algorithm.

PAGE 52

41 vi) at block interfaces the displacement is interpolated as the average between neighboring blocks values, in order to maintain the point-to-point interface. The main advantage of this algorithm is its ability to handle large rotational and translational deformations, while maintaining a good grid quality near the body; as one can see in Fig. 2-4(A). A) B) Figure 2-4. Grid quality for the grid re-arrangement algorithm. A) Close to the body. B) Overall domain. The disadvantages of this algorithm are the large skewness near the outer boundary (Fig. 2-4(B)) imposed by the non-movement condition and the approximation of dx, dy, d quantities, leading to poor grid quality after a long time. This problem is resolved for periodic movements by replacing the deformed grid with the original one after each cycle. 2.5.3 Grid Re-arrangement Algorithm Based on Spring Analogy Concept To overcome the main disadvantage of the grid re-arrangement algorithm described in section 2.5.2 and of the master-slave algorithm (section 2.5.1), a new algorithm based on the spring-analogy concept was developed by Tang et al. (2006). This algorithm propagates the perturbation due to body movement along a grid line as described in

PAGE 53

42 section 2.5.1, but also takes into account the angular deformation weighted with the original grid orthogonality information. This way, the cross-over encountered in the original master-slave algorithm (Fig. 2-2) is eliminated. Table 2-1. Moving boundaries techniques. Advantages and disadvantages. No. Grid movement algorithm Advantages Disadvantages Observations 1. Master-slave (sec. 2.5.1) Can handle large translational displacements Displacement can be controlled by varying the stiffness Cannot handle large rotation deformations Used for MAV studies small vertical disp. 2. Weight function (sec 2.5.2) Can handle large translational and rotational displacements. Displacement can be controlled by varying the stiffness parameter (p) Error accumulation after long time Used for 2-D flapping wing studies (large rotations). 3. Spring analogy (sec. 2.5.3) Can handle large translational and rotational displacements No current generalization for arbitrary multiblock grids. Used for 2-D flapping wing studies (large rotations). Extensive computations also proved that the grid quality is maintained in time. The main advantages and disadvantages of the grid movement routines are summarized in Table 2-1. 2.6 Membrane Structural Solver The main characteristic of membrane structures is their large deformations implying a nonlinear stress-strain relation. Green and Adkins (1960) were the first to derive a general nonlinear theory for membranes. The membrane theory is based on two main assumptions (Jenkins and Leonard, 1991): the membrane tension stiffness is much larger than its bending stiffness and the ratio between the membrane thickness and the radius of the curvature is small (therefore the strain-displacement is decoupled from the curvature tensor). Two approaches exist to model the response of a membrane structure. The first approach idealizes a plate or shell structure as a membrane away from the boundaries and neglects the stress couples in the interior region. This model can be

PAGE 54

43 described by Foppl-von Karaman theory (Jenkins and Leonard, 1991; Jenkins, 1996). The second approach assumes that the structure cannot sustain stress couples anywhere. This approach is used in the current work and was also adopted by Lian et al. (2002) in their simulation of a membrane wing. The studies of a two-dimensional membrane with fixed leading and trailing edges were originally conducted by Thwaites (1961) and Nelsen (1963). In these studies, an inextensible membrane surrounded by steady, two-dimensional, irrotational flow was considered. Under these assumptions, and for small camber and incidence angle, the problem can be linearized and expressed in an integral form known as the sail equation (Newman, 1987). A comprehensive review of earlier works related to membrane wing aerodynamics can be found in Newman (1987). Unlike the two-dimensional case, the three-dimensional membrane model introduces several complicating factors. In the two-dimensional sail equation, the tension is a single constant. In three-dimensional membranes, the tension is defined as a biaxial tension along the lines of principal stress (Jackson and Christie, 1987). Also, the geometric and material properties may vary along the spanwise direction and need to be described in detail. Another complicating factor is that the membrane can be compressed and wrinkled when one of the principal tensions vanishes. A finite element analysis of the static equilibrium of an inflated membrane undergoing large deformations was presented by Oden and Sato (1967). A review of the earlier work on the dynamic analysis of membranes can be found in Jenkins and Leonard (1991). An update of the state of the art models in membrane dynamics is presented by Jenkins (1996). Verron et al. (2001) have studied, both numerically and experimentally,

PAGE 55

44 the dynamic inflation of a rubber-like membrane. Ding et al. (2003) numerically studied partially wrinkled membranes. The three-dimensional membrane model used in this work to study the membrane dynamics has been proposed by Lian et al. (2002). The model gives good results for membrane dynamics with large deformations but cannot handle the wrinkle phenomenon that occurs when the membrane is compressed (Lian, 2003a). The membrane material considered obeys the hyperelastic Mooney-Rivilin model (Mooney, 1940). A brief review of their membrane model is given next. 2.6.1 Mooney-Rivlin Model The Mooney-Rivlin model is one of the most frequently employed hyperelastic models because of its mathematical simplicity and relatively good accuracy for reasonably large strains (less than 150%) (Mooney, 1940). For an initially isotropic membrane, a strain energy function, W, can be defined as (Green and Adkins, 1960): 321,,IIIWW (2-12) where I 1 I 2 and I 3 are the first, second, and third invariants of the Green deformation tensor, C. For an incompressible material, when I 3 = 1, the strain energy is a function of I 1 and I 2 only, and a linear expression can be written for the membrane strain energy: 332211IcIcW (2-13) where c 1 and c 2 are two material constants. A material that obeys equation (2-13) is known as a Mooney-Rivlin material. For an initially isotropic, incompressible membrane, the general stress-strain relation is written as:

PAGE 56

45 CICS212112cIccp (2-14) where S is the second Piola-Kirchhoff stress tensor, p is the hydrostatic pressure, and I is the 3 identity matrix. The membrane stress field is essentially assumed to be two-dimensional and therefore, the stress normal to the deformed membrane surface is negligible with respect to the stress in the tangent plane (Oden and Sato, 1967). Under this assumption, the deformation matrix C(t) and the stress matrix S(t) can be written as: 0000)()(0)()()(;)(000)()(0)()()(221212113322121211tStStStSttCtCtCtCtCtSC (2-15) The hydrostatic pressure is determined by the condition that S 33 = 0, and the formula is: 232321212cIccp (2-16) where 3 is defined by 0333)()(hthtC (2-17) in which h(t) and h 0 are the membrane thickness in the deformed and undeformed configurations, respectively. 2.6.2 Solution of the Dynamic Equations Using the discretized form of the principle of virtual work in the Lagrangian description, the expression for the Green-Lagrange strain tensor, and the second Piola-Kirchhoff stress tensor, the dynamic equation of motion for one triangular element in global coordinates can be written as: exteeeeTeeFFDMDint (2-18)

PAGE 57

46 where e is the virtual work for one element, D e is the virtual nodal displacement vector, M e is the mass matrix, D is the nodal acceleration vector, F is the internal force vector due to membrane deformation, and F is the external load vector. Subscript e indicates values for one triangular element. The expressions for the element mass matrix, M e inte exte e the internal force vector, F, and the external nodal load, can be found in Lian (2003a). Equation (2-18) gives the virtual work for a single element. If we sum all elements in global coordinates and take into account that the virtual work of the whole system is zero, the governing equation for the membrane structure dynamic response under external loads is: inte exteF exttFFDMint)( (2-19) where M is the positive definite mass matrix (which remains constant) D(t) is the nodal displacement vector for all nodes, D is the nodal acceleration vector, F )(t int is the internal elastic resisting force of the membrane due to its deformation, and F ext is the nodal external load vector. To solve the dynamic membrane equation (equation (2-19)), either an implicit or an explicit scheme can be used for marching in time. The explicit scheme is computationally economical for each time step and requires less storage. However, the time step needs to be small in order to satisfy the stability requirements. The implicit scheme requires more computational time and storage per time step because the system matrix is iteratively updated at each time step to account for the nonlinear effects. But, generally, they are unconditionally stable, and there is no restriction in the choice of time step size. In this work, the implicit Wilsonmethod (Wilson, 1973) is used. It is a one-step method which

PAGE 58

47 is easier to code and behaves well in fluid-structure interaction problems. The Wilsonmethod has a dissipative character; therefore a smaller time step than that required by the overall accuracy criterion is used to avoid excessive dissipation. A larger time step tends to cause more numerical dumping as shown by Lian (2003a). Also, the method is unconditionally stable for linear problems and is second-order accurate in time. For details and illustrations of the membrane model and its computations results, see Lian et al. (2002), and Lian (2003a). 2.7 Coupling between Fluid Solver and Structural Solver During flight, the surrounding flow exerts aerodynamic forces on the membrane structure which leads to deformations in the wing shape. At the same time, the change in shape of the solid surface modifies the flow around the wing. To capture the interaction between fluid and structure, the governing equations for both fluid and structure need to be solved at the same time. However, structure and fluid formulations have different mathematical and numerical properties, and depending on the problem complexity, different techniques need to be developed to provide coupling between the computational fluid dynamics (CFD) and computational structure dynamic (CSD) methods. One approach is to treat the flow equations and the structural equations as one system of equations (Bendiksen, 1991) and solve the system of equations simultaneously on a single grid. However, the presence of multiple time scales in the fluid and structure system makes a unified treatment inefficient in handling the computational stiffness. To adequately handle different characteristics of the flow and structure interaction, many researchers choose to solve the fluid and structure equations separately and couple them through a synchronization process that involves subiterations (Smith and Shyy, 1997; Gordnier and Fithen, 2003; Zwaam and Pranata, 2002; Lian and Shyy, 2005; Lian

PAGE 59

48 et al., 2002). The fluid solver and the structural solver run independently on their own computational grid with their own choice of time step. Then a subiteration between the two solvers is employed at each time step to couple the structure and fluid system of equations and avoid phase lag. Under the moving boundary framework for fluid-structure interaction problems, the CFD grid needs to be automatically re-generated after the shape change. Since different grids are used by the structure and fluid solvers, methods to transfer the information between the CFD and CSD grids need to be employed (Smith and Hodges, 2000). Lian et al. (2002), and Lian and Shyy (2005), in their study of the membrane wing dynamics, have adopted a thin plate spline (TPS) interpolation method to exchange information between the fluid and structure solver. The same method is also used in this work. The thin plate spline interpolation is a global interpolation method which is invariant with respect to rotation and translation, and thus is a good choice for interpolation on moving surfaces. The interpolation function for the one-dimensional case can be written as: NiiiixxxxxH12log)( (2-20) where H(x) is the interpolated displacement, i is the undetermined coefficient, N is the number of structural nodes on the interface, and x i is their location. Once the coefficients i are found based on the structural grid locations, the nodal displacement vector defined on the CSD grid, D, is mapped to the CFD grid through an interpolation matrix, G, GDX (2-21) where X is the corresponding displacement on the CFD grid. Similarly, the pressure forces can be mapped from the CFD grid to the CSD grid.

PAGE 60

49 Start the time loo p t = 0 Fluid Solver Evaluate pressure pi by solving the Navier-Stokes e q uations on the CFD g rid. Interface Transfer the new pressure pi to the CSD grid throu g h TPS inter p olation Structural solver Based on interpolated p evaluate the nodal dis p lacement vector D i Interface Map the displacements from CSD to CFD g rid usin g TPS inter p olation: X i =GD i Re-meshing Based on the interpolated displacements u p date the CFD g rid Re-compute the CFD grid metric terms Update the Jacobian to satisfy the GCL N ext time ste p t = t + t STOP Read the CFD and CSD grids Subiteration to cou p le fluid and structure solvers i = Figure 2-5. Generic algorithm for fluid-structure interaction In order to conduct time-accurate computations for fluid-membrane dynamics, iterations need to be done at each time instant to enforce the fluid and structure coupling. A generic algorithm for fluid-structure interaction problems with separate fluid and

PAGE 61

50 structure solvers that exchange information between the CFD and CSD grid through interpolation, and uses subiterations to enforce the fluid-structure coupling, is shown in Fig. 2-5.

PAGE 62

CHAPTER 3 FIXED WING MICRO AIR VEHICLES 3.1 Overview Micro air vehicles with a maximum dimension of 15 cm or less and a flight speed of 10-20 m/s are of interest to both military and civilian applications. Equipped with a video camera or a sensor, these vehicles can perform surveillance and reconnaissance, targeting, and bio-chemical sensing at a remote or otherwise hazardous location (Shyy et al., 1999). With the rapid progress made in structural and material technologies, and miniaturization of power plants, communication, visualization, and control devices, several groups have developed successful MAVs (Pornsin-Sirirak et al., 2000; Grasmeyer and Keennon, 2001; Ifju et al., 2002; Jones and Platzer, 2003; Jones et al., 2004). With its reduced size, the MAV gains certain favorable scaling characteristics, including low inertia and reduced stalling speed (defined as the minimum speed at which the wing can produce sufficient lift force for flight) (Shyy et al., 1999). However, in terms of aerodynamics, control, range and maneuverability, many outstanding issues exist. MAVs fly in a low Reynolds number regime (10 4 ~10 5 ) due to their low flight speed and small dimensions. This is often accompanied by laminar boundary layer separation, transition, and reattachment, which leads to a decrease in the aerodynamic performance (Shyy et al., 1999). Valuable information in these regards is reviewed by Lissaman (1983) and Tani (1964), as well as addressed in several proceedings by Mueller (1998, 2000). Active control approaches with micro-electronic-mechanical systems (MEMS) for MAVs are discussed by Gad-el-Hak (2001). 51

PAGE 63

52 In the development of MAVs, there are three main approaches: (1) using flapping wings for generating lift and thrust; (2) employing rotating wings; (3) using fixed wings for generating lift. From the design point of view, fixed wing configuration is simple in concept and easy to implement. On the other hand, flapping wings, which show superiority over fixed wings in terms of lift generation and maneuverability at low Reynolds number, is unfortunately difficult to put into practice due to complex wing kinematics and power requirements. The simplicity in design of the fixed wing MAV, as a miniature of large airplane wings, is counterbalanced by a deterioration in performance because of the presence of separation bubbles, as its operating Reynolds number drops to the range of 10 4 to 10 5 Also, the efficiency of conventional propellers diminishes with decreasing diameter. Flapping wing design, which is predominantly selected as the optimal approach by nature, compensates the complex kinematics by simultaneously providing both lift and propulsion. Figure 3-1(A) demonstrates a fixed wing design by Ifju and co-workers at the University of Florida (Ifju et al., 2002). This design features a unique flexible wing, which is advantageous to the MAV flight because it can adapt to wind gusts, offering a smoother platform. Also, experimental studies show that the flexible wing has higher stall angles than its counterpart rigid wing. The membrane wing consists of a leading edge spar and six chord-wise battens, which are made of unidirectional carbon fiber prepreg laminates. The membrane material is bonded to the spar and battens. While the membrane is flexible, and cannot sustain bending moments, the carbon fiber is rigid, can sustain bending moments and provides structural support.

PAGE 64

53 The other approach, inspired by nature, uses flapping wings to produce both lift and thrust. This concept has been practiced in the designs of Pornsin-Sirirak et al. (2000) and is shown in Figure 3-1(B). A mixed design is presented by Jones et al. (2004). Their MAV has a fixed wing used to generate the lift force, while using small flapping wings to generate thrust. A B Figure 3-1. Fixed and flapping wing MAVs. A) University of Florida 15 cm (6) MAV design with a flexible wing. B) CalTech/AeroVironment Microbat flapping wing MAV. In this chapter we focus on the aerodynamics of both rigid and flexible fixed wing MAVs, based on the designs of Ifjus group at the University of Florida. The steady, rigid wing aerodynamics is presented first, followed by the unsteady flow over a flexible wing.

PAGE 65

54 3.2 Rigid Wing Micro Air Vehicles 3.2.1 Wing Geometry Figure 3-2 shows the wing shapes designed by Ifju and co-workers for their MAVs and are used as a base for the present numerical simulations. Figure 3-2(A) shows an older design of a single curvature airfoil wing, which has a 15 cm (6) span, a root chord of 13.3 cm (5.24) and an area of 160 cm 2 (24.8in 2 ). A more recent design is shown in Figure 3-2(B). The wing has a span of 12.7 cm (5), a root chord of 10.4 cm (4.09) and an area of 110.4 cm 2 (17.11in 2 ). The wing has two distinct regions: a central panel with almost constant airfoil profile and no dihedral angle; and an outer panel with variable cross-section and a dihedral angle of almost 7 o The double curvature airfoil cross-section offers better transversal (pitch) stability characteristics for flying wing airfoils. A) B) Figure 3-2. MAVs wing shapes. A) 15cm (6) MAV wing shape. B) 12.7 cm (5) MAV wing shape. 3.2.2 Computational Setup To compute the aerodynamic performance of the MAV wing, the steady state, Navier-Stokes equations for the incompressible flow are solved as described in Chapter 2. A pressure-based algorithm is used to solve the three-dimensional Navier-Stokes equations written in curvilinear coordinates. Both first-order and second order upwind

PAGE 66

55 schemes are employed for convection terms and second-order central difference schemes are adopted for pressure and viscous terms. For the flow over the rigid MAV wings, the physical problem under consideration is the steady state flow over a wing in an unbounded domain. For all simulations, a structured, multi-block grid is employed in the computations. The outer boundaries are placed far away from the wing to reduce the boundary influence on the flow near the wing. The schematic of the computational domains for the 6 wing and 5 wing is presented in Figure 3-3. INLET OUTLET INLET OUTLET (-6c, 5c, -6c) (9c, 5c, -6c) (-6c, -5c, 6c) A B (9c, -5c, 6c) Figure 3-3. Computational domain for the MAV wing. A) 15 cm (6) wing. B) 12.7 cm (5) wing. In this work and in similar studies (Lian and Shyy, 2005; Lian et al., 2002; Lian et al., 2003c), only the wing is taken into consideration, and no effects from the fuselage and propeller are accounted for. Since the freestream velocity is along the chord-wise direction and no propeller is modeled, only half of the wing is modeled (based on the symmetry condition). A uniform velocity is specified at the inlet boundaries, while for the outlet boundaries a zero-gradient condition is imposed. On the wing surface, the noslip condition is specified.

PAGE 67

56 3.2.3 Tip Vortex Effect on the Rigid Wing Aerodynamics On a finite wing, the pressure difference between the high-pressure region below the wing surface and the low-pressure region above the wing generates a circulatory motion around the wing tip that affects the wing aerodynamics. The total drag coefficient for a finite wing at subsonic speed can be written as (Anderson, 1989): A ReCCCCLFDPDD2,, (3-1) where C D,P is the drag coefficient due to pressure, C D,F is the drag coefficient due to skin friction, e is the span efficiency factor (which is less than one), AR is the aspect ratio, and A ReCCLiD2, is the induced drag coefficient due to the existence of tip vortices. Equation (3.1) demonstrates that the induced drag varies as the square of the lift coefficient; at high angles of attack the induced drag can constitute a substantial portion of the total drag. The MAV wings presented by the UF MAV group (Ifju et al., 2001; Ifju et al., 2002; Albertani et al., 2004), seen in Figure 3-2, have a low aspect ratio of approximately 1.4; the wing shape is dictated by the need to maximize the wing area, and hence the lift, for a given dimension. Equation (3.1) illustrates that as AR is decreased, the induced drag increases. Therefore, for these low-aspect wings, it is important to investigate the tip vortex effects on the wing aerodynamics. In general, tip vortex effects are two-folded: Tip vortex causes downwash that decreases the effective angle of attack and increases the drag force (Anderson, 1989). Tip vortex forms a low-pressure region on the top surface of the wing, which provides additional lift force (Mueller and DeLaurier, 2003).

PAGE 68

57 Lian et al. (2003c), and Lian and Shyy (2005) performed numerous computations to study the incompressible, laminar flow over the rigid 15 cm (6) MAV wing (Fig. 3-2(A)). The relative position of the wing with respect to the outer-boundaries of the computational domain is presented in Figure 3-3(A). Based on the root chord length and a free-stream uniform velocity of 10m/s, the Reynolds number is 9x10 4 They studied the flow over the wing at angles of attack of 6 o 15 o 27 o 39 o and 51 o Figure 3-4 shows the vorticity contours and pressure on the wing surface as they evolve from 6 o to 27 o angle of attack. At 6 o angle of attack, the tip vortices are clearly visible, but their strength is modest and affects only a small portion of the wing. As the angle of attack increases to 27 o the tip vortices are stronger and the low-pressure area on the wings upper surface becomes larger. Figure 3-4. Streamlines, vorticity and pressure on the 6 MAV wings upper surface as a function of the angle of attack. A) 6 o angle of attack. B) 27 o angle of attack. Adopted from Lian et al. (2003c). At higher angle of attack (45 o ), the flow separates from the leading edge, and the low-pressure region associated with the tip vortex core increases, which helps to maintain the increase in lift force.

PAGE 69

58 Lian et al. (2003c) and Lian and Shyy (2005) concluded that for low-aspect ratio wings, the tip vortex has an important effect on the wing aerodynamics. On one hand, the tip vortex increases the lift at high angle of attack, similar to the vortex lift phenomenon observed on delta wings (Torres and Mueller, 2001). On the other hand, the tip vortex increases the induced drag, and therefore decreases the aerodynamic performance. If we look at the MAV as a surveillance platform, it is desirable to have good flight efficiency, which translates to larger values of the lift-to-drag ratio. To accomplish this objective we try to maintain/increase the lift while decreasing the induced-drag. Various methods to reduce the induced drag by decreasing the tip vortex effects are described in the literature and confirmed by actual applications to aircraft wing design (LaRoche and Paffy, 1996). In this study, the focus is on the implication of placing endplates at the wing tip, as this is the simplest method from a manufacturing point of view. In this work, we study the wing tip vortex effect on the rigid 15 cm (6) MAV wing (Fig. 3.2(A)) and the influence of an endplate on the overall flow. Since a curved endplate deteriorates the overall aerodynamic performances (Viieru et al., 2003), a modified wing platform that allows straight endplates was created. Three wing geometries have been studied: (1) the original wing (Fig. 3-5(A)) with dimensions specified in section 3.2.1; (2) the modified wing (Fig. 3-5(B)) with a span of 14 cm (5.5) and a wing area of 155 cm 2 (24 in 2 ), while the root chord has the same length as the original wing; (3) the modified wing with endplates attached. The endplate has a length of 4.4 cm (1.73) and a height of 3.4 cm (1.34). The modified wing with endplates is shown in Fig. 3-5(C).

PAGE 70

59 Original wing shape B C A Figure 3-5. Wing shape geometry. A) Baseline wing. B) Modified wing. C) Endplates location on the modified wing. The arrow indicates the flow direction. The computational domain is identical to that presented in Fig. 3-3(A). The boundary conditions are described in section 3.2.2, with the addition of noslip boundary conditions assigned to the endplate surfaces. The free-stream velocity is 10 m/s for all cases, which gives a Reynolds number of 9x10 4 based on the root chord length. Lian and Shyy (2005) have performed computations for the same original wing geometry (Fig. 3-5(A)) using three different grid sizes ranging from 180,000 nodes (41x31 on the wing surface) at the coarse level to 2.3 million nodes (101x61 on the wing surface) at the fine level. They concluded that the coarse grid solution under-predicts the lift-to-drag ratio (C L /C D ) by 0.5% and form drag force by less than 0.7%, compared to the finest grid solution. A similar result was also obtained for the lift force. However the skin friction drag is sensitive to grid density, and the computation based on the coarse mesh overestimates it by 10%. Since in the presently considered parameter range, the form drag is at least one order of magnitude higher than the viscous drag, we use the coarse grid in our numerical simulation to save computational time. In Fig. 3-6, the vorticity magnitude contours behind the wing are plotted in planes perpendicular to the x-axis, i.e., the stream-wise direction. The streamlines are also displayed to clearly show the wing tip vortex. It can be seen that both wings without the

PAGE 71

60 endplate show a strong wing tip vortex, while for the modified wing with endplates, the tip vortex core is smaller and dissipates faster. B A C Figure 3-6. Vorticity magnitude contours behind the wing at 6 o angle of attack. A) Original wing. B) Modified wing. C) Modified wing with endplates. Arrow indicates the flow direction. The endplate basically deters the flow from the high-pressure region on the lower wing surface from reaching the low-pressure region on the upper wing surface. In the absence of the endplate, the fluid is accelerated towards the wing tip by the vortex cores low-pressure zone present on the upper wing surface, as indicated by the streamlines plotted in Fig. 3-7(A), (B). The low-pressure zone associated with the tip vortex core can be clearly identified in the plot of the span-wise pressure variation on the wings upper surface (Fig. 3-9(A) and Fig. 3-10(A)).

PAGE 72

61 A B C Figure 3-7. Pressure contours on the wings upper surface and streamlines slightly above the wing at 6 o angle of attack. A) Original wing. B) Modified wing. C) Modified wing with endplates. Pressure is in [N/m 2 ]. The bending of the streamlines towards the wing tip in the absence of the endplates, caused by the pressure difference between the high-pressure zone (below the wing surface) and the low-pressure zone (above the wing surface) can also be observed by looking at the streamlines slightly below the wing surface, as shown in Fig. 3-8(A) and (B). One may also note that the tip vortex is stronger in the case of the modified wing without endplates than in the case of the original wing. This points to a more significant pressure drop on the lower wing surface, and can be observed in Fig. 3-8(B). The decrease of the high-pressure area due to the circulatory motion in the absence of the

PAGE 73

62 endplate can also be recognized by looking at the span-wise pressure coefficient on the lower wing surface, which is plotted in Fig. 3-9(B) and Fig. 3-10(B). A B C Figure 3-8. Pressure contours on the wings lower surface and streamlines slightly below the wing at 6 o angle of attack. A) Original wing. B) Modified wing. C) Modified wing with endplates. Pressure is in [N/m 2 ]. The endplate decelerates the flow near the wing tip. On one hand, this decrease in velocity reduces the pressure drop on the upper wing surface corresponding to the vortex core (Fig. 3-9(A) and Fig. 3-9(B)), reducing the vortex lift. On the other hand, a lower velocity slightly below the wing increases the high-pressure area since more momentum is transferred to the wing as pressure, instead of being shed as vorticity at the wing tip. The increase in the high-pressure zone on the lower wing surface in the presence of the endplate can be clearly seen from the span-wise pressure coefficient on the lower wing

PAGE 74

63 surface plot (Fig. 3-9(B) and Fig. 3-10(B)), and also from the pressure contours on the lower wing surface shown in Fig. 3-8. 0 0.1 0.2 0.3 0.4 0.5 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 Z / Local span-Cp Original wing, NO endplatesModified wing, NO endplatesModified wing, WITH endplates 0 0.1 0.2 0.3 0.4 0.5 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 Z / Local span-Cp Original wing, NO endplatesModified wing, NO endplatesModified wing, WITH endplates A) B) Figure 3-9. Pressure coefficient on the wing surface at x/c = 0.34 and 6 o angle of attack. A) Upper wing surface. B) Lower wing surface. 0 0.1 0.2 0.3 0.4 0.5 -0.3 -0.25 -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 Z / Local span-Cp Original wing, NO endplatesModified wing, NO endplatesModified wing, WITH endplates 0 0.1 0.2 0.3 0.4 0.5 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 Z / Local span-Cp Original wing, NO endplatesModified wing, NO endplatesModified wing, WITH endplates A) B) Figure 3-10. Pressure coefficient on the wing surface at x/c = 0.53 and 6 o angle of attack. A) Upper wing surface. B) Lower wing surface. To get a better picture of the endplate effects, the pressure difference between lower and upper wing surfaces was integrated along the local chord direction at different span-wise locations. The result is the lift and drag estimation in the span-wise distribution. Figure 3-11 shows the span-wise lift distribution. It can be noticed that the endplate, by blocking the flow around the tip, reduces the downwash. Therefore the

PAGE 75

64 effective angle of attack is increased, leading to an increase in lift towards the wing tip. The figure shows that when the endplates are attached, the lift is increased when compared with the wing shape without the endplate. 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0. 5 1.5 2 2.5 3 3.5 4 4.5 Z / Original Wing SpanLift / Span [N/m] Original wing, NO endplatesModified wing, NO endplatesModified wing, WITH endplates Figure 3-11. Span-wise lift distribution at 6 o angle of attack. The modified wing with endplates produces almost the same amount of lift as the original wing. However, if we look at the span-wise drag distribution plotted in Fig. 3-12, the modified wing with or without endplate offers the lowest drag for almost 75% of the half wingspan (starting from the root). After computing the total lift and drag by integrating the pressure difference on the wing surface, we conclude that for an angle of attack of 6 o an improvement in lift-to-drag ratio of 10.1% for the modified wing configuration with endplates is observed as compared to the baseline configuration (original wing, no endplate) (Viieru et al., 2004). This improvement is based mainly on the reduction in drag offered by the modified wing shape, since the lift force profiles are similar.

PAGE 76

65 For 15 o angle of attack, the modified wing with endplates shows an increase of 1.4% in lift-to-drag ratio compared to the baseline configuration. This is achieved by an important reduction of drag, while lift is decreasing with by 5%. 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0. 5 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 Z / Original Wing SpanDrag / Span [N/m] Original wing, NO endplatesModified wing, NO endplatesModified wing, WITH endplates Figure 3-12. Span-wise drag distribution at 6 o angle of attack. 3.3 Flexible Wing Micro Air Vehicles As mentioned before, flexible wings increase the aerodynamic performance of fixed and flapping wings MAVs, but due to the complexities involved in the fluid-structure problem, less rigorous work has been considered in this regard. Waszak et al. (2001) performed an experimental study of a flexible wing MAV. Lian et al. (2003c) numerically simulated the flow over a flexible wing. They showed that a MAV with a flexible wing yields better performance than that of a rigid wing because of its shape adaptation and delay to stall. Ho et al. (2003) experimentally demonstrated that flexible wings have better lift production for flapping flight. However, no decisive conclusion has been drawn on how the wing flexibility enhances the aerodynamic performance. The influence of the wing flexibility of a 12.7 cm (5) MAV at 6 o angle of attack is studied next. The geometry is identical with that presented in section 3.2.1 (Fig. 3-2(B)).

PAGE 77

66 The central wing panel is rigid, while the outer panels are covered with a flexible material that has a uniform thickness of 0.2 mm and a density of 1200 kg/m 3 (Ifju et al., 2001; Ifju et al., 2002; Lian et al., 2002). The flexible part has a rigid leading edge spar and two flexible carbon-fiber battens, as shown in Fig. 3-13(B). The battens are modeled (in the structural solver) as membrane material with increased density to simulate the dynamic response. For the finite element structural solver (chapter 2), an unstructured mesh with 598 nodes and 1076 elements was created. A coarse CFD grid with 210,000 nodes was used, since the overall aerodynamic performances are not very sensitive to grid size (Lian and Shyy, 2005) and significant computational time can be saved. Both CFD and FEM grids are presented in Fig. 3-13. Figure 3-13. Surface grid for the 5 MAV wing. A) CFD coarse mesh. B) FEM mesh for the flexible part of the wing indicating the boundary conditions for the structural solver. Based on a free-stream uniform velocity of 10 m/s and the root chord length, the Reynolds number for this case is 7.1x10 4 For the unsteady computations, a time step of 2x10 -4 s is used to ensure that membrane vibration time scales are captured. Two sub

PAGE 78

67 iterations were used to ensure the time synchronization as shown by Lian et al. (2003c). The fluid-structure framework is presented in chapter 2. Even though the flow stream is uniform, the membrane exhibits self-initiated vibrations, as one can see in Fig. 3-14. This was also observed by Lian et al. (2003c) in their computations of the flow over a 15 cm (6) MAV with flexible wing. 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 -0.1 -0.05 0 0.05 0.1 0.15 0.2 Vertical dispacement 5in MAVtime [sec]Vert disp [cm] near wingtipbetween the battens Figure 3-14. History of the vertical displacements of two points, one near the wing tip and the other one between the two battens. It can be clearly seen that the point near the tip exhibits larger displacements than the point between the two battens. The probable explanation of this fact is the presence of the strong circulatory motion of the tip vortex. Also, the battens rigidity restrains the movement of the points between them. The vertical displacement contours are visualized in Fig. 3-15 at different time instants. The first snapshot shows how the membrane starts to inflate in the region of largest pressure difference. The ripples in the membrane, especially between the battens, can be observed in the next time snapshot. This is correlated with the periodic behavior seen in Fig. 3-14.

PAGE 79

68 Figure 3-15. Vertical displacement contours at various times. The displacement is in [cm]. To assess the overall aerodynamic performances, the pressure and shear stress are integrated over the wing surface at each time step. In Fig. 3-16 the lift and drag coefficients history are presented for the 12.7 cm (5) MAV wing at 6 o angle of attack. 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 Lift coefficient history, 5in MAV, AoA = 6 degtime [sec]CL 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0.04 0.045 0.05 0.055 0.06 0.065 0.07 0.075 0.08 Drag coefficient history, 5in MAV, AoA = 6 degtime [sec]CDA) B) Figure 3-16. Aerodynamic forces history for the 5 flexible wing MAV. A) Lift coefficient. B) Drag coefficient.

PAGE 80

69 Next, the time averaged aerodynamic performance is compared with the steady-state, rigid wing at the same angle of attack and with experimental values measured in the wind tunnel for the same rigid wing geometry (Viieru et al., 2004). The results are compared in Table 3-1. Table 3-1. Aerodynamic performances for flexible and rigid wing. Parameter Computational Flexible wing (time average) Computational Rigid wing Experimental Rigid wing Lift coefficient (C L ) 0.37 0.40 0.42 Drag coefficient (C D ) 0.064 0.076 0.085 Lift / Drag ratio (C L /C D ) 5.78 5.26 5.06 In the computation of the lift and drag coefficients, both pressure and viscous forces are taken into consideration. The flexible wings averaged lift coefficient is smaller, but comparable with the rigid wings computational and experimental values. This indicates that for small angles of attack, the aerodynamic performances are not very different. In terms of drag coefficient, it was found to be smaller for the flexible wing, leading to a better lift-to-drag ratio. The best explanation for the drag reduction for the flexible wing is the adaptation to the flow, resulting most probably in a reduction of the tip vortex effects and therefore a reduction in the induced drag. The small difference between the aerodynamic performance of the rigid and flexible wing was also noticed by Lian et al. (2003c) in his computations for the 6 MAV wing. However, experimental work has proven that the advantages of flexible wings come into play at high angles of attack. Generally, membrane wings show a higher stall angle than rigid wings (Waszak et al., 2001; Albertani et al., 2004). Due to complicated flow behavior at high angles of attack for this range of Reynolds numbers, mainly flow separation and transition from

PAGE 81

70 laminar to turbulent boundary layer, numerical investigation of the behavior of the flexible wings at high angles of attack remains an open topic. 3.4 Summary and Conclusions for Fixed Wing Aerodynamics In this chapter, we numerically investigated the flow over rigid and flexible fixed wings. The steady-state, laminar, incompressible Navier-Stokes equations were solved using a pressure based solver to study the flow over 6 and 5 wingspan wings at Reynolds numbers ranging from 7x10 4 to 9x10 4 at different angles of attack. For the flexible wing computations, the unsteady Navier-Stokes solver was coupled with a finite element structural solver to investigate the fluid-structure interaction problem of a flexible wing in a uniform flow. First, the effect of the tip vortex on the overall aerodynamics of a low aspect ratio 15 cm (6) MAV wing at 6 o and 15 o angles of attack was studied. To reduce the induced drag generated by the tip vortex, the endplate effect was also studied for both the same wing configuration and for a modified wing that allows for the installation of a straight endplate. The following observations were made: The modified wing shape reduces the drag, with or without endplate. However the lift is also reduced due to a stronger vortex and a reduced area. The endplate parallel to the freestream velocity increases lift without a significant increase in drag. The increase in lift is closer to the tip, where the endplate reduces the flow travelling from the lower surface to the upper surface of the wing, reducing the downwash, and increasing the effective angle of attack. Consequently the lift increases. The effectiveness of the endplate diminishes as the angle of attack increases, due to a stronger wing tip vortex. Secondly, the impact of wing flexibility was studied on a 12.7 cm (5) wingspan wing at 6 o angle of attack and Reynolds number of 7.1x10 4 The results showed that the

PAGE 82

71 membrane wing vibrates, even in a steady free-stream. For small to moderate angles of attack, and before stall limit, the time-averaged lift and drag coefficients of the membrane wing seem to be comparable with those of a rigid wing. However, the flow behavior at high angles of attack for this range of Reynolds numbers, which departs form the pure laminar assumption, and from the high Reynolds number turbulent models, presents increased difficulties for the numerical investigation of the flexible wings and remains an open topic for research.

PAGE 83

CHAPTER 4 AERODYNAMICS OF LOW REYNOLDS NUMBER HOVERING AIRFOILS Even though three-dimensional effects are critical for MAVs, experiments and computations in two-dimensions continue to provide important insight into flapping wing aerodynamics. From the computational point of view, 2-D simulations require less computational power, and grid sensitivity and convergence studies can be performed more readily (compared to 3-D cases). Wang (2000a, 2000b) demonstrated that the aerodynamics of 2-D wings could help explain some of the enhanced lift observed in insect flight. In this chapter, the flow structures around a hovering 2-D elliptic airfoil are studied. First, the viscous force computation is validated for the Reynolds number range encountered in insect flight. Second, the results for a hovering airfoil are compared with similar numerical and experimental work to validate the current fluid solver. Finally, the lift generation mechanisms are compared for two hovering modes, and the effect of the Reynolds number on the aerodynamics is investigated. 4.1 Validation of Viscous Force Computation Low Reynolds number flows are characterized by large viscous forces that become closer in magnitude to the pressure forces, and therefore have a large contribution to the overall lift and drag forces. Generally, the shear stress at the wall is computed as the product of the velocity gradient in the normal direction at the wall surface and the viscosity: wallwallun (4-1) 72

PAGE 84

73 In the present approach, equation 4-1 is discretized using a first order scheme: tantan P wallwalluuntan() (4-2) where, P u)wall is the tangential velocity at the center of the first cell from the wall, is the tangential velocity of the wall (non-zero for moving boundary problems) and n is the normal distance from the wall to the center of the first cell near the wing surface. tan(u To validate the present formulation, two cases are studied. First, the steady flow over a flat plate at zero angle of attack at different Reynolds numbers representative of insect flapping flight is computed. The evaluated friction coefficient is compared with analytical results. Secondly, the flow over an oscillating flat plate in a quiescent medium is simulated in order to assess the viscous force computation for moving walls. 4.1.1 Steady Flow over a Flat Plate at Zero Angle of Attack To simulate the steady flow over a thin plate of chord c = 1.0, the Navier-Stokes equations are solved on a grid with 100 points along the streamwise and 60 points along the cross-stream direction. The distance from the wall to the first cell center is 5x10 -4 c, which guarantees a sufficient number of points in the boundary layer for high Reynolds number flows. The numerical simulations are performed for Reynolds numbers (based on chord length and free-stream velocity), of 10, 10 2 10 3 and 10 4 The second order upwind scheme was employed for convection terms, while second order central difference schemes are adopted for the pressure and viscous terms. A no-slip boundary condition is imposed on the airfoil surface.

PAGE 85

74 In Fig. 4-1, the numerical and analytical velocity distribution in the boundary layer for a flat plate at different Reynolds numbers is plotted. The analytical velocity distribution is defined as (Schlichting, 1979): ()uUf (4-3) where U is the free-stream velocity, f() is the solution of Blasiuss ordinary differential equation (Blasius, 1908), and /()yUx is a dimensionless coordinate. 0 1 2 3 4 5 6 7 8 9 0 0.2 0.4 0.6 0.8 1 1.2 U velocity profile, Re = 10,000u/U x/c = 0.99x/c = 0.50x/c = 0.20Blasius 0 1 2 3 4 5 6 7 8 9 0 0.2 0.4 0.6 0.8 1 1.2 U velocity profile, Re = 1,000u/U x/c = 0.99x/c = 0.50x/c = 0.20Blasius 0 1 2 3 4 5 6 7 8 9 0 0.2 0.4 0.6 0.8 1 1.2 U velocity profile, Re = 100u/U x/c = 0.99x/c = 0.50x/c = 0.20Blasius A) B) C) Figure 4-1. Numerical and analytical velocity distribution in the boundary layer along different locations on a flat plate and different Reynolds numbers. A) Re=10 4 B) Re=10 3 C) Re=10 2 Here /()yUx is the dimensionless vertical coordinate.

PAGE 86

75 Figure 4-1 shows good agreement between the numerical results and the analytical formulation, especially for high Reynolds numbers. As the Reynolds number decreases, the computed velocity distribution moves further away from the analytical distribution. This is expected, since the Blasiuss equation is based on Prandtls (1904) boundary layer equations, which implied the high Reynolds number flow assumption. The analysis of the flow over a flat plate using Blasius equation is restricted to semi-infinite plate, because of the parabolic nature of the Prandtls boundary layer equations that cannot account for upstream changes in shear stress initiated by the trailing-edge of a finite-length plate. The discrepancy between the analytical and numerical velocity distribution is visible at locations near the trailing-edge of the plate, as one can observe in Fig. 4-1 for all Reynolds numbers. Figure 4-2 also shows that as the Reynolds number decreases, the friction coefficient given by Blasius departs from the numerical solution. To take into account second order effects, generalizations of Prandtls boundary layer equations were developed. For the case of the flat plate, Stewartson (1974) and Messiter (1970) found an expression that improves the prediction of the skin-friction coefficient for low Reynolds number. The analytical expressions for skin-friction coefficient, defined for 2-D cases as are presented in Table 4-1. 2/(0.5)fCDragcU Table 4-1. Numerical and analytical friction coefficient values for flat plate. Case Reynolds number C f Blasius (1904) 1.328/ReflC C f Messiter (1970) 7/81.328/Re2.668/ReflC l C f present numerical results 1 104 0.0133 0.0141 0.0142 2 103 0.0420 0.0483 0.0486 3 102 0.133 0.180 0.182 4 10 0.420 0.776 0.785

PAGE 87

76 Figure 4-2 shows a very good agreement between the improved analytical solution of Messiter (1970) and the numerical results for a large range of Reynolds numbers, therefore validating the viscous force computation method employed in the solver. 101 102 103 104 0.01 0.02 0.04 0.06 0.1 0.2 0.4 0.6 1 Friction coefficient vs Reynolds numberReCf BlasiusMessiterStream Figure 4-2. Numerical and analytical friction coefficient C f for a flat plate versus Reynolds number. 4.1.2 Flow over an Oscillating Flat Plate in a Quiescent Environment In the previous sub-section, the viscous force computation is validated for steady flows at Reynolds numbers from 10 to 10 4 Next, the concept needs to be validated for the case of moving walls (wall velocity is non-zero). If we consider an infinite plate oscillating along the x-axis, the wall velocity is defined as: 0(0,)cos(uytUt ) (4-4) where U 0 is the maximum velocity and 2 f ; f being the oscillation frequency. The integration of equation 4-4 gives the plate displacement: 0()sin();/aahthtwithhU (4-5)

PAGE 88

77 The analytical velocity field for this motion is given by Stokes (1851): 0/exp()cos()uUt (4-6) where /2y is a non-dimensional cross-stream coordinate. Following the definition given in equation 4-1, the analytical wall shear stress for an oscillating plate is given by (White, 1991): 0sin(/4)wallUt (4-7) The unsteady flow over an oscillating plate of length 1 is solved on a grid with 50 points along the chord and 80 points in the vertical direction using the Navier-Stokes solver. The distance from the wall to the first cell center is 10 -3 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 Wall shear stress for oscillating platetime/T, T= w, uw/10 w computeduwall/10 computedw analytic soluwall/10 analytic sol Figure 4-3. Skin-friction and wall velocity for an oscillating plate. Oscillating frequency /21/f amplitude h a =0.5 and Reynolds number =10 3 The wall velocity is scaled to an order of magnitude comparable to the wall shear stress. Continuous lines = present computation, Symbols = analytical solution for infinite plate.

PAGE 89

78 For the case studied, 2 results in a frequency f of 1/ and a period T of To ensure a sufficient number of points for one oscillation cycle, a time step size of 0.01 is chosen. Based on the maximum velocity U 0 and chord length, the Reynolds number is 10 3 In Fig. 4-3, the numerically predicted wall shear stress on the finite plate is compared with analytical values for an infinite plate given by equation 4-7. The figure shows a good agreement between numerical and analytical values. The small phase shift between numerical and analytical results can be explained by the effects of the leading and trailing edges of a finite plate. In Fig. 4-3 the scaled wall velocity is also plotted to contrast the lag between maximum shear and maximum velocity. t/T = 0.45 t/T = 0.29 t/T = 0.21 t/T = 0.05 Figure 4-4. Velocity distribution above an oscillating plate. Oscillating frequency /21/f amplitude h a =0.5 and Reynolds number =10 3 Continuous lines = present computational results. Symbols = analytical solution for infinite plate.

PAGE 90

79 The theoretical velocity distribution above the plate, given by equation 4-6, is plotted, along with the computed velocity field, in Fig. 4-4 for different time instants during one half-period. The skin friction and velocity profiles presented in Fig. 4-3 and Fig. 4-4 show a good agreement between theoretical and computed results. Good agreement between the analytical solution for steady flow over a flat plate, as well as for an oscillating plate is attained, validating the viscous force computation methodology employed in the fluid solver for stationary or moving walls. 4.2 Kinematics of 2-D Airfoil in Hovering Flight In hovering flight, insects and small birds draw the still ambient fluid and move it downwards in a succession of vortex dipoles (or rings in 3-D) that carry the fluid momentum and generate the lift necessary to stay aloft. The hovering flight plays an important role for insects in feeding. Different hovering patterns are visualized in nature. Figure 4-5 shows the figure pattern observed in hovering hummingbirds. A B Figure 4-5. Hummingbird in hovering flight. A) Forward stroke, B) Backward stroke. Picture taken by Wei Shyy. In an extensive survey, Weis-Fogh (1973) noticed that most insects, such as fruit flies, bees, and beetles, employ symmetric back-and-forth strokes near a horizontal stroke plane, as seen in Fig. 4-6(A). Weis-Fogh referred to this style as normal hovering. Some other insects, like true hoverflies and dragonflies employ asymmetric strokes along an inclined plane. In this case the lift force is primarily generated during the downstroke

PAGE 91

80 (Wang, 2005). The dragonflys inclined plane hovering schematic is presented in Fig. 4-6(B). Since the normal hovering is employed by most of the insects, it receives most of the focus in experimental and computational investigations (Dickinson et al., 1999; Ellington et al., 1996; Sun and Tang, 2002a, 2002b; Wang et al., 2004a). In this study, we focus on two hovering modes, the water treading and the normal hovering modes. The kinematics of these hovering modes are presented in the next section. B A Figure 4-6. The schematics of two hovering styles. A) Hovering using a horizontal plane. B) Hovering using an inclined stroke plane. Sketched from Wang (2005). 4.2.1 Kinematics of Water Treading Mode The water treading hovering mode is one of the cases that Liu and Kawachi (1998) used to validate their finite volume algorithm. The hovering mode studied is based on the so-called water treading mode as defined by Freymuth in his experiments (Freymuth, 1990) and consists of a translational component (plunging), given by: ()sin(2)ahthkt (4-8) and a rotational component (pitching), given by: 0()sin(2)atk t (4-9) In equation 4-8, h a is the plunging amplitude and k is the reduced frequency of the sinusoidal oscillation. In equation 4-9, 0 is the initial angle of attack, a is the pitching

PAGE 92

81 amplitude, and is the phase difference between the plunging and pitching motion. The reduced frequency k is defined by: 2/2/2refakfcUc h (4-10) where c is the airfoil chord length, f is the oscillation frequency, and U 2refa fh is the reference velocity (equal to the maximum plunging velocity). It should be noticed that the reduced frequency, by definition, varies with the inverse of the stroke amplitude. By setting the reference velocity (U ref ) to one, the relationship between the reduced frequency, k, and the frequency, f, is: f=k/c, and the stroke period is T=c/k. The Reynolds number is defined as Re /refUc where is the kinematic viscosity of the fluid. A schematic of the airfoil movement is shown in Fig. 4-7. ha-ha-a X Y ha -ha+a X Y Forward stroke Backward stroke Figure 4-7. Schematics of the water treading hovering mode. The specific parameters for the airfoil movement will be specified for each individual case in this study that employs water treading hovering mode.

PAGE 93

82 4.2.2 Kinematics of Normal Hovering Mode In this work a normal hovering mode based on a sinusoidal motion (Wang et al., 2004) is also studied. The equations describing the elliptic airfoil motion are: ()sin(2)ahthkt t (translation) (4-11) 0()sin(2)atk (rotation) (4-12) where h a is the stroke amplitude, 0 is the initial angle of attack, a is the pitching amplitude, is the phase difference between the plunging and pitching motion, and k is the reduced frequency as defined in equation 4-10. A schematic of the motion is shown in Fig. 4-8. -a ha -ha Y X Forward stroke (t) +a ha -ha Y X Backward stroke (t) Figure 4-8. Schematics of normal hovering mode. Based on the above equations, the parameters describing the airfoil motion will be specified for each case studied in this work. 4.3 Water Treading Hovering Mode Liu and Kawachi (1998) placed the airfoil in a body-centered inertial frame of reference that undergoes the hovering motion about a fixed axis at origin. The plunging

PAGE 94

83 motion is implemented in the boundary conditions at the open boundary stencils. The O-type grid used in their work has different sizes, varying from 71x71 to 145x99 grid points, in order to investigate the effect of the grid resolution. The minimum grid spacing adjacent to the solid wall is 2.4x10 -3 c. In this work, an O-type grid with different resolutions was also employed. At the coarse level, the grid has 90x70 nodes, while at the fine level the grid has 250x120 nodes. The computational domain for the coarse grid is shown in Fig. 4-9. The minimum grid spacing adjacent to the airfoil varies from 2.4x10 -3 c (coarsest grid) to 1.0x10 -3 c (finer grid). Figure 4-9. The coarse multi-block grid used in the current work. The grid has 90 points around the airfoil and 70 points in the radial direction. The distance from the airfoil to the outer boundary is 15 chords. The unsteady, laminar, incompressible, Navier-Stokes equations are solved using the pressure based algorithm described in chapter 2. To handle the large airfoil displacements and rotations, a moving grid algorithm is coupled with the flow solver.

PAGE 95

84 4.3.1 Grid and Time Increment Sensitivity Analysis The spatial accuracy of the present algorithm is examined by employing three grid size levels. A refinement of the time step size was also performed on the intermediate grid (named Dgr1). The hovering mode is described by equations 4-8 and 4-9 with the same kinematics parameters as in the water treading hovering mode presented earlier. Table 4-2 provides a summary of the cases run, along with grid and time increment details. A schematic of the movement is presented in Fig. 4-7. Table 4-2. Grid and time increment details for the water treading hovering cases. Here y represents the distance from the wall to the first grid point and t is the time step size. Re=1,700, h a =1.0c, 0 =0 o a =42 o k=0.5, = -90 o Case Grid t Computational time (hrs) (single processor) Number of time steps 1. Dgr0 (90x70, y=2.4x10 -3 c) 0.01 7.5 (Intel Itanium 1.3GHz) 5600 (~ 8 periods) 2. Dgr1 (190x80, y =1.0x10 -3 c) 0.01 31.4 (Intel Itanium 1.3GHz) 5600 (~ 8 periods) 3. Dgr2 (250x120, y =1.0x10 -3 c) 0.01 64.3 (Intel Itanium 1.3GHz) 5600 (~ 8 periods) 4. Dgr1 (190x80, y =1.0x10 -3 c) 0.02 13.3 (Dec Alpha 850MHz) 2800 (~ 8 periods) 5. Dgr1 (190x80, y =1.0x10 -3 c) 0.005 44.8 (Intel Itanium 1.3GHz) 11200 (~8 periods) The force coefficients are obtained by normalizing the forces by: 20.5()Vt where is the density and 2()0.5refU 2 Vt (Liu and Kawachi, 1998). Figure 4-10 and 4-11 show the lift and drag coefficient history for three different grid resolutions. The results show a good agreement between solutions and indicate that reasonable grid independence was obtained. The lift and drag forces attain periodic behavior after almost 5 periods. The highest lift peak values are observed for the finest grid (250x120). To study the effect of time step size with a good spatial accuracy, the intermediate grid, named Dgr1 (190x80), is used. Figures 4-12 and 4-13 show the lift and drag coefficient history for time increments varying from 0.005 to 0.02

PAGE 96

85 5 5.5 6 6.5 7 7.5 8 -2 -1 0 1 2 3 4 5 time/TCLLift coefficient history for t=0.01 Dgr0 ( 90x70)Dgr1 (190x80)Dgr2 (250x120) Forward stroke Backward stroke T=2 Figure 4-10. Lift coefficient for three periods and different grid sizes and t=0.01. h a =1.0c, 0 =0 o a =42 o k=0.5 and Re=1,700. 5 5.5 6 6.5 7 7.5 8 -4 -3 -2 -1 0 1 2 3 4 time/T, T=2 CDDrag coeffcient history for t=0.01 Dgr0 ( 90x70)Dgr1 (190x80)Dgr2 (250x120) Forward stroke Backward stroke Figure 4-11. Drag coefficient for three periods and different grid sizes and t=0.01. h a =1.0c, 0 =0 o a =42 o k=0.5 and Re=1,700.

PAGE 97

86 5 5.5 6 7 7.5 8 -2 -1 0 1 2 3 4 5 time/TCLLift coeffcient history on grid Dgr1 (190x80) t=0.02t=0.01t=0.005 Forward stroke Backward stroke 6.5 T=2 Figure 4-12. Lift coefficient for three periods and different time increments. h a =1.0c, 0 =0 o a =42 o k=0.5 and Re=1,700. 5 5.5 6 6.5 7 7.5 8 -4 -3 -2 -1 0 1 2 3 4 time/T, T=2 CDDrag coeffcient history on grid Dgr1 (190x80) t=0.02t=0.01t=0.005 Forward stroke Backward stroke Figure 4-13. Drag coefficient for three periods and different time increments. h a =1.0c, 0 =0 o a =42 o k=0.5 and Re=1,700.

PAGE 98

87 As observed in Fig. 4-10, the lift coefficient reaches a plateau close to the maximum pitch angle (the middle of forward/backward stroke), and then a sharp increase in lift is noticed as the airfoil starts pitching down. This is captured well for time increments of 0.005 and 0.01; a smaller time step better captures the vortex shedding phenomenon responsible for the sudden decrease in lift. Using a time step size of 0.02, the sharp increase in lift following the pitch-down of the airfoil is not captured very well. Based on the results for different spatial and temporal resolutions, we conclude that the combination of an intermediate grid (190x80) and a time step size of 0.01 gives satisfactory balance between numerical accuracy, the flow structures, and computational cost. The results also indicate that reasonable grid independence was obtained, taking into account the unsteadiness of the flapping flight. 4.3.2 Comparison with Similar Work Based on the grid and time step size refinement study, the results obtained on the intermediate grid (190x80) with a time increment of 0.01 are compared with the work of Liu and Kawachi (1998). Figure 4-14(A) shows the lift and drag coefficients history for one period as obtained by Liu and Kawachi, while Fig. 4-14(B) shows the current computed values. The lift and drag coefficients obtained in the present numerical simulation correspond well to those of Liu and Kawachi, and further confirm the capabilities of the Navier-Stokes solver used in this work. 4.3.3 Lift Generation Mechanism for Water Treading Hovering Mode. In this section, the lift generation mechanism is investigated for this particular hovering mode. The flow analysis is based on the last computed period in order to avoid the transient effects and to allow the force variation to reach a periodic pattern. The data

PAGE 99

88 obtained using the finest grid (250x120) and a time increment of 1x10 -2 is used for this analysis. Eight time instants are selected for analysis during one hovering period and are identified in Fig. 4-15. The airfoil position and pitching angle, as well as the translational and rotational velocities, are plotted in Fig. 4-15(A). Fig. 4-15(B) shows the lift and drag coefficient variations. 19 20 21 22 23 24 25 -4 -3 -2 -1 0 1 2 3 4 Lift and Drag coefficients for one periodTimeLift & Drag coefficients CDCL B) A) Figure 4-14. Lift and Drag coefficients for one period. h a =1.0c, 0 =0 o a =42 o k=0.5 and Re=1,700. A) Liu and Kawachi (1998) (figure obtained via personal communication). B) Present computational results. At the beginning of the forward stroke, the airfoil accelerates from rest and starts to pitch-up rapidly. Fig.4-15(B) shows that most of the lift and drag force is generated when the airfoil travels from time instant A to B. Between time instant B and C, the airfoil travels close to the maximum translation velocity and at an almost constant maximum pitching angle (Fig. 4-15(A)). Figure 4-15(B) shows that, during this period, the lift curve reaches a plateau. The clockwise vortex shed earlier (Fig. 4-16(A) at time A) accelerates the flow over the upper surface of the airfoil, and creates a low-pressure zone (Fig. 4-16(A) and (B), time B).

PAGE 100

89 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 -3 -2 -1 0 1 2 3 4 Lift and Drag cots for one periodtime/CL, CD CDCL A B C D F G H E 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 Kinematics for "water treading" hovering modeti2 h(t), (t), dh/dt, d /dt alpha [rad]stroke displacementtranslational velocityrotational velocity F o r wa r d st r o k e Backward stroke A) B) efficien me/T, T= T, T=2 Figure 4-15. Kinematics of water treading hovering mode and the aerodynamic forces for one complete stroke. h a =1.0c, 0 =0 o a =42 o k=0.5 and Re=1,700. A) kinematics. B) lift and drag coefficients. The non-dimensional time instants marked on the figure correspond to: A=0.03, B=0.19, C=0.31, D=0.46, E=0.55, F=0.67, G=0.79 and H=0.95. After reaching the maximum translational velocity and pitching angle at time/T = 0.25 (see Fig. 4-15(A)), the airfoil starts decelerating and pitching down, and a counter-clockwise leading edge vortex starts to form (Fig. 4-16(A), at time C). The pressure drop inside the vortex core increases (Fig. 4-16(B) at time C): therefore the lift reaches the maximum peak during the half-stroke (~ time/T = 0.3). Once the leading edge vortex is

PAGE 101

90 shed (Fig. 4-16(A) at time D), the lift decreases sharply and becomes negative (Fig. 4-15(B)), close to the end of the forward stroke. A B C D E F G H A B C D E F G HA) B) Figure 4-16. Flow field snapshots of the water treading hovering mode. h a =1.0c, 0 =0 o a =42 o k=0.5 and Re=1,700. A) Vorticity contours; Red=counter-clockwise vortices, Blue = clockwise vortices. B) Pressure contours; Red = high pressure, Blue = low pressure. The capital letters on each snapshot correspond to the time instants defined in Fig. 4-12. The symmetry of the lift and drag variation during the forward and backward stroke observed in Fig. 4-15(B) is also supported by the symmetric flow structures (Fig. 4-16(A) and (B), at times F, G and H). The maximum lift peak during the backward stroke is seen around time/T = 0.78 (see Fig. 4-15(B)), after the airfoil starts decelerating and rotates down, similar to what was observed in the forward stroke. The symmetric movement of

PAGE 102

91 the airfoil is reflected in the symmetry of the force variation, and suggests that no major interaction between the airfoil and the vortex structures generated in the previous periods is present. The transient effects due to the start-up condition of the hovering motion are attenuated quickly, as the force variation becomes periodic after 3-4 cycles. 4.4 Normal Hovering Mode 4.4.1 Grid Refinement Study The normal hovering mode, where the wing moves in a figure pattern, is the hovering mode most employed by insects and small birds, and is the focus of the following section. For the present 2-D investigation, the movement of the airfoil undergoing normal hovering is described by equation 4-11 (translation) and 4-12 (rotation). The schematics of the normal hovering mode are presented in Fig. 4-8. The unsteady, laminar, incompressible, Navier-Stokes equations are solved in an O-type grid domain. The spatial accuracy of the present algorithm is examined by employing three grid sizes. At the coarsest level, the grid has 90 points around the airfoil and 70 points in the radial direction. The computational domain for the coarse grid is shown in Fig. 4-9. Table 4-3 shows the summary of cases run, including the grid and time increment details. Table 4-3. Grid and time increment details for the normal hovering cases. Here y represents the distance from the wall to the first grid point and t is the time step size. Re=100, h a =1.4c, 0 =90 o a =45 o k=1/2.8, = -90 o Case Grid t Computational time (hrs) (single processor) Number of time steps 1. Dgr0 (90x70, y=2.4x10 -3 c) 0.01 11.26 (Dec Alpha 850MHz) 7000 (~ 7 periods) 2. Dgr1 (190x80, y =1.0x10 -3 c) 0.01 103.68 (Dec Alpha 850MHz) 7000 (~ 7 periods) 3. Dgr2 (250x120, y =1.0x10 -3 c) 0.01 379.28 (Dec Alpha 850MHz) 7000 (~ 7 periods) For all cases mentioned in Table 4-3, a second order upwind scheme was employed for convection terms, while second order central difference scheme are adopted for the pressure and viscous terms.

PAGE 103

92 The forces are normalized by the maxima of their corresponding quasi-steady forces, as described in Wangs work (Wang et al., 2004). For 2-D computations, the quasi-steady translational lift force is found to be well approximated by:, where 20.5qsLqsLuC 1.2sin(2)LqsC 2max0.5maxLu (Wang et al., 2004). Therefore, the lift force is normalized by: qsLqsC In Fig. 4-17, the lift coefficient history for the last three periods is plotted. At a Reynolds number of 100, a periodic pattern is noticed after 4 periods. The coarse grid (named Dgr0) result features asymmetric lift peaks during half-strokes for one period, as shown in Fig. 4-17. The results obtained on the intermediate (named Dgr1) and fine grid (named Dgr2) show no clear asymmetry in lift peaks during the half-strokes. The very small difference observed between fine and intermediate grids indicate that a grid independent solution was obtained. Forward stroke Backward stroke 4 4.5 5 5.5 6 6.5 7 -1 -0.5 0 0.5 1 1.5 2 Lift coefficient for normal hoveringtime/T, T=2.8 CL Dgr0 ( 90x70)Dgr1 (190x80)Dgr2 (250x120) Forward stroke Backward stroke Figure 4-17. Lift coefficient of the normal hovering mode for three periods and different grid sizes and t=0.01. h a =1.4c, a =45 o k=1/2.8, Re=100.

PAGE 104

93 To conclude, the grid refinement study shows that the intermediate grid (190x80) gives sufficient resolution for an airfoil undergoing normal hovering at a Reynolds number close to 100. 4.4.2 Comparison with Similar Computational and Experimental Results In this section, the current computational results are qualitatively compared with the experimental and computational results of Wang et al. (2004), to observe the similarities in the flow structures and aerodynamic force pattern. The flow data was obtained using a dynamically scaled robot fly, described in Dickinson and Gtz (1993), Dickinson et al. (1999), and Birch and Dickinson (2003). The wing is cut to the planform of a Drosophila wing and is made from 2.25 mm thickness Plexiglas. When attached to the force sensor arm, the total wing length is 0.25 m. The wing and arm apparatus is lowered in a Plexiglas tank (1m x 1m x 2m) filled with mineral oil (density = 0.88x103 kg m -3 kinematic viscosity = 115cSt). To measure the flow structures, digital particle image velocimetry (DPIV) was used. The light sheets were parallel with the wing chord and positioned at 0.65R, where R is the wing span. More details about the experimental setup can be found in Wang et al. (2004). For the present computational setup, the translation motion is described by equation (4-11), and the rotational motion is given by equation (4-12). The schematics of the airfoil motion in the hovering mode are shown in Fig. 4-8. The kinematics parameters for both the experimental and computational setup are specified in Table 4-4. Based on the grid sensitivity analysis, the intermediate grid (190x80) is used to generate the present results with a good spatial accuracy.

PAGE 105

94 The lift force presented in Fig. 4-18 is normalized with the maximum quasi-steady force as described in Wang et al. (2004) and section 4.4.1. Table 4-4. Kinematics parameters for the experimental and computational setup. Here c is the chord length and the Reynolds number is based on the maximum velocity and the chord length. Case Plunging amplitude (h a ) Initial pitching angle ( 0 ) Pitching amplitude ( a ) Phase difference () Frequency Reynolds number 1. Experimental (Wang et al., 2004) 1.4c 90 o 45 o -90 o f=0.25Hz 75 2. Computational (Wang et al., 2004) 1.4c 90 o 45 o -90 o f=0.25Hz 75 3. Computational (Present) 1.4c 90 o 45 o -90 o k=1/2.8 75 4. Experimental (Wang et al., 2004) 2.4c 90 o 45 o -90 o f=0.25Hz 115 5. Computational (Present) 2.4c 90 o 45 o -90 o k=1/4.8 115 Figure 4-18 shows a similar computed force pattern and lift coefficient values with the computational and experimental measurements of Wang et al. (2004). Computational results of Wang et al. (2004) (Fig. 4-18(A)) show asymmetric force peaks during one stroke, which was not clearly observed in the present computations (Fig. 4-18(B)) and the experimental measurements. However, one must be cautious when comparing 2-D and 3-D lift coefficients. The experimental measurements indicate less asymmetric time histories between forward and backward strokes than Wang et al.s computational results. The present computational results exhibit similar trends. However, the maximum and minimum values are somewhat different. Furthermore, in Fig. 4-19, the computed vorticity field is compared with the experimental flow visualization data of Wang et al. (2004). The non-dimensional stroke amplitude employed in the experiment and present computational setup is h a /c = 2.4 (case 4 and 5 in Table 4-4) and the Reynolds number is

PAGE 106

95 115. The computed (Fig. 4-19(A) and (C)), and experimentally measured (Fig. 4-19(B) and (D)) vorticity fields show a good agreement in terms of predicted flow structures. 0 0.5 1 1.5 2 2.5 3 3.5 4 -1 -0.5 0 0.5 1 1.5 2 Lift coefficient historytime/TCL CL B) A) Figure 4-18. Lift coefficient history for the normal hovering mode. A) Results of Wang et al. (2004) for symmetric rotation. h a =1.4c, and Re=75. Red = experimental measurements. Blue = computational results. Reprinted with permission from The Company of Biologists Ltd. B) Present computational results for h a = 1.4c, a =45 o ,k=1/2.8 and Re=75. The normal hovering study demonstrated the capability of the present numerical framework to predict, with a certain degree of confidence, the aerodynamic forces generated on a 2-D elliptic airfoil in hovering motion. It is concluded that a grid with 190 points around the airfoil and 80 in the radial direction offers a good spatial resolution and a reasonable computational time for the Reynolds number selected and, therefore, this grid will be used in future cases.

PAGE 107

96 t/T=6.9 t/T=6.8 t/T=6.7 t/T=6.6 t/T=6.5 t/T=6.4 t/T=6.3 t/T=6.2 t/T=6.1 t/T=6.0 B) D) A) C) Figure 4-19. Vorticity plot for normal hovering mode with h a /c=2.4 and Re=115. Red = counter-clockwise vortices, Blue = clockwise vortices. For all cases, the stroke reversal is symmetric. (A, C) Computed vorticity for the sixth stroke. (B, D) DPIV data in a 2D slice at 0.65R. The frames are measured during the fourth stroke. The color scales for computed and measured vorticity plots do not correspond to the same contour values and should be viewed as a qualitative comparison. Experimental vorticity contours from Wang et al. (2004). Reprinted with permission from The Company of Biologists Ltd. 4.5 Lift Generation Mechanisms in Normal and Water Treading Hovering Modes In this work, the aerodynamic force generation by the same 15% thickness elliptic airfoil undergoing two different hovering modes is studied. First, we consider the water treading mode (Freymuth, 1990) described in section 4.2.1 and governed by equations 4-8 (translation) and 4-9 (rotation). Next, we consider the normal hovering mode

PAGE 108

97 presented in section 4.2.2, described by the equations in 4-10 and 4-11. To compare the two hovering modes, the kinematics parameters are selected to be consistent with each other, and are presented in Table 4-5. Table 4-5. Kinematic parameters for water treading and normal hovering modes. The Reynolds number for both cases is 100. Hovering mode Initial angle of attack 0 Pitching amplitude a Stroke amplitude h a /c Reduced frequency k Phase difference water treading hovering 0 o 45 o 1.4 1/2.8 -90 o normal hovering 90 o 45 o 1.4 1/2.8 -90 o Figure 4-7 shows the schematics of the water treading hovering mode. The schematic of the normal hovering mode is presented in Fig. 4-8. For both cases, the intermediate grid (190x80) was used and a time step size of 0.01 is chosen for the unsteady computations. The Reynolds number for this case is 100 based, on the maximum plunging velocity and the airfoils chord. Figure 4-20 shows the lift and drag coefficients during one complete stroke for water treading and normal hovering modes. The aerodynamic forces in the figures above are normalized by 0.5, where U 2refU ref = 1. In the case of the water treading hovering mode, for the first half of the forward stroke the airfoil accelerates and pitches-up. During this interval, the lift increases constantly (Fig. 4-20(A), 0 to c), and the unsteady dynamics prevent the flow separation from the upper side of the airfoil, even at high angles of attack, as indicated by the vorticity contours plotted in Fig. 4-21a to c. The maximum lift peak is obtained close to the middle of the half-stroke, when the pitch angle reaches the maximum value (Fig, 4-20c). After this stage, the airfoil starts to decelerate, and pitches-down. The flow separates and a large re-circulation bubble forms

PAGE 109

98 on the upper side of the airfoil (Fig. 4-21d and e), leading to a decrease in lift to the minimum value (Fig. 4-20(A), at time e). The same pattern is repeated for the backward stroke. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 Lift coefficient for two hovering modes at Re=100time/TCL "normal" hovering mode"water treading" hovering mode a b c d e f g h 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 Drag coefficient for two hovering modes at Re=100time/TCD "Normal" hovering mode"Water treading" hovering mode a b c d e f g h A) B) Figure 4-20. One cycle force history for two hovering modes. h a =1.4c, a =45 o k=1/2.8 and Re=100. A) Lift coefficient. B) Drag coefficient. The selected normalized time instants are: a) 0.08, b) 0.17, c) 0.25, d) 0.31, e) 0.45, f) 0.60, g) 0.80, h) 0.94.

PAGE 110

99 d normal hovering c b a h normal hovering g f e water treading water treading Figure 4-21. Vorticity contours for two hovering modes. h a =1.4c, a =45 o k=1/2.8 and Re=100. Red = counter-clockwise vortices, Blue = clockwise vortices. The flow snapshots (a to h) correspond to the time instants defined in Fig. 4-20. At the beginning of the forward stroke of the normal hovering mode, the airfoil starts accelerating and pitches-down. The rotation of the airfoil accelerates the flow around the leading and trailing edges, creating a suction zone on the upper side of the airfoil, while the high-pressure stagnation area on the lower side is increased due to fluid driven from the surroundings by the previously formed vortex (Fig. 4-21a). This combination of low and high-pressure areas leads to an increase in lift at the beginning of the stroke (Fig. 4-20(A), at time a). As the airfoil rotates downward and accelerates, the fluid is accelerated towards the trailing edge and the high-pressure stagnation area decreases (Fig. 4-21b). The lift also decreases, reaching a local minimum at time/T ~0.17

PAGE 111

100 for the forward stroke and 0.57 for the backward stroke, as shown in Fig. 4-20(A). At the middle of the half-stroke, the airfoil travels at almost constant pitching angle. A re-circulation bubble attached to the airfoil forms on the upper surface (Fig. 4-21c, d, g) and helps increase the lift and drag to their maximum values (time/T ~ 0.3 and 0.8) during one complete stroke (Fig. 4-20(A) and (B), at d and g). After the maximum pitching angle and translation velocity are obtained (time/T= 0.25 and 0.75) during one half-stroke, the airfoil decelerates and pitches up, leading to flow separation on the upper side of the airfoil (Fig. 4-21e and h). The detachment of the large vertical structure from the upper airfoil surface, combined with rapid deceleration, decreases the circulation. Therefore, the lift coefficient drops to its minimum value (Fig. 4-20(A) at times e and h). The force coefficient history for water treading and normal hovering modes indicates differences in the lift generation mechanism. For both hovering modes, the lift force reaches its maximum value when the airfoil moves near maximum velocity at almost constant maximum pitching angle. Similar maximum lift peak values (Fig. 4-20(A) at times d and g) and flow structures (Fig. 4-21d and g) are observed in this particular time interval (mid-stroke), suggesting a similar lift generation mechanism. The vorticity contours (Fig. 4-21) indicate that the flow is either attached, or have a small re-circulation bubble on the upper side of the airfoil. Therefore, the delayed-stall mechanism is mainly responsible for generating most of the lift force. The delayed-stall mechanism was also observed during wing translation by Sun and Tang (2002) in their simulation of a Drosophila wing in hovering flight for a similar Reynolds number range. While the delayed-stall is the main lift generation mechanism in the case of the water treading hovering mode, for the normal hovering mode the local lift peaks at

PAGE 112

101 the beginning of the half-strokes point out that a wake-capturing mechanism is also a contributing factor. Hence, depending on the detailed kinematics, the lift generation mechanisms at a Reynolds number of 100 exhibit different physical mechanisms. The averaged lift coefficient for both cases is computed as the summation of the lift coefficient over the last three periods divided by the total time. For the water treading hovering mode an average lift coefficient of 0.77 is obtained, while for the normal hovering mode the average lift coefficient is 0.56, suggesting that the water treading mode performs better at Re=100 with the given kinematic parameters. The significant role of viscosity at low Reynolds numbers reduces the interaction between vortex structures generated during previous strokes, as reflected by the almost symmetric maximum peaks for lift and drag as seen in Fig. 4-20(A) and (B). To investigate the Reynolds number effect on the aerodynamic forces and the flow structure, we have computed the hovering aerodynamics of the water treading mode at Re=100 and Re=1,700. Based on the same kinematics of the Re=100 case, the aerodynamics of the water treading mode is assessed. The kinematics and flow parameters for these cases are summarized in Table 4-6, and the airfoils motion schematic is presented in Fig. 4-7. Table 4-6. Parameters for water treading hovering mode employed for Reynolds number effect study. Hovering mode Initial angle of attack 0 Pitching amplitude a Stroke amplitude h a /c Reduced frequency k Phase difference Reynolds number Re 1. Water treading 0 o 45 o 1.4 1/2.8 -90 o 100 2. Water treading 0 o 45 o 1.4 1/2.8 -90 o 1,700 The lift coefficients for the water treading hovering mode at Reynolds 100 and 1,700 are plotted in Fig. 4-22. The lift coefficient peaks are noticeably higher for the

PAGE 113

102 Reynolds number of 1,700 than that for a Reynolds number of 100. While the force patterns between the two Reynolds numbers are similar, the higher Reynolds number case exhibits less symmetry between forward and backward strokes. The vorticity (Fig. 4-23(A) and (C)) and pressure distribution on the airfoil surface (Fig. 4-24(A)) indicate a symmetric flow field for a Reynolds number of 100. +a 4 4.5 5 5.5 6 6.5 7 -0.5 0 0.5 1 1.5 2 2.5 Lift coefficient for "water treading" mode at Re=100 and 1,700time/TCL Re=100Re=1,700 a b c d Backward stroke Forward stroke ha -ha Figure 4-22. Lift coefficient for the water treading mode. h a = 1.4c, a =45 o k=1/2.8 and a Reynolds number of 100 and 1,700. The selected normalized time instants are: a) 6.25, b) 6.48, c) 6.77, d) 6.97. The pressure distribution on the airfoil surface, plotted in Fig. 4-24, show that near the maximum lift peaks, the high pressure stagnation area on the lower side of the airfoil is almost similar in both shape and magnitude for the two Reynolds number studied. However, on the upper side of the airfoil, the mild variation of the pressure gradient for

PAGE 114

103 the low Reynolds number case (Fig. 4-24(A), time a and c) suggests that the flow is attached, while for the high Reynolds number (Fig. 4-24(B), time a and c) the low-pressure area near the leading edge indicates a re-circulation zone corresponding to the formation of the leading edge vortex. This low-pressure area is responsible for most of the high lift peak values seen in the case of a Reynolds number of 1,700 (Fig. 4-22, times a and c). In the lower Reynolds number regime, the viscous effect is capable of offsetting the asymmetry effects caused by the start-up condition of the flapping motion. At the higher Reynolds number regime, the convective effect dominates. Hence the initial condition persists to form the asymmetric patterns between the forward and backward stroke generated by the interaction of the vortex structures history and vortex structures that are dissipated slowly for high Reynolds numbers cases. a b c d D) C) B) A) Figure 4-23. Vorticity contours for the water treading mode. h a = 1.4c, a =45 o k=1/2.8. Red = counter-clockwise vortices, Blue = clockwise vortices. (A, C) Reynolds number = 100. (B, D) Reynolds number = 1,700. The flow snapshots (a to d) correspond to the time instants defined in Fig. 4-22.

PAGE 115

104 c c a chord chord chord a chord B) A) Figure 4-24. Pressure distribution on the airfoil surface for the water treading mode. h a = 1.4c, a =45 o k=1/2.8. A) Reynolds number = 100. B) Reynolds number = 1,700. The flow snapshots (a, c) correspond to the time instants defined in Fig. 4-22. 4.6 Summary and Conclusion for the Aerodynamics of Hovering Airfoils In this chapter, the flow over an elliptic airfoil in hovering motion under different flow parameters was numerically investigated. The unsteady, laminar, incompressible Navier-Stokes equations were solved using a structured, pressure-based solver to study the aerodynamics of hovering airfoils at Reynolds numbers ranging from 100 to 1,700. The flow solver was coupled with a grid re-arrangement algorithm (Tang et al., 2006) to create a numerical framework capable of accommodating large translational and rotational displacements, characteristic of flapping flight. First, the water treading (Freymuth, 1990) hovering mode is investigated. The numerical results, for stroke amplitude of 1 chord and pitch amplitude of 42 o at a

PAGE 116

105 Reynolds number of 1,700, show good agreement with the numerical results of Liu and Kawachi (1998), validating the present approach. The flow analysis indicates that the delayed-stall mechanism is responsible for most of the lift generation during one complete stroke. Second, the aerodynamics of the normal hovering mode, as defined by Weis-Fogh (1973), at low Reynolds numbers is investigated. The present numerical results are compared with the experimental and numerical results of Wang et al. (2004) for stroke amplitudes of 1.4 and 2.4 chords at Reynolds numbers of 75 and 115. The computational results of Wang et al. (2004) show clear asymmetric lift force patterns during one stroke. The present results exhibit similar trends, but the asymmetry of the lift force peaks during one complete stroke is less noticeable, similar to the experimental measurements (Wang et al., 2004). To get more insight in the lift generation mechanisms, the aerodynamics of airfoils undergoing hovering flight at a Reynolds number of 100, based on the water treading and normal hovering modes is investigated. The delayed-stall mechanism is found to be responsible for generating the maximum lift peaks for both hovering modes. However, for the normal hovering mode, a wake-capturing mechanism is identified at the beginning of each stroke, consistent with the observations of Wang et al. (2004) for similar kinematics. Finally, the effect of the Reynolds number on the aerodynamics of the water treading-based hovering is studied. The delayed-stall is the main lift generation mechanism at both Reynolds numbers studied. The less symmetric force patterns at the higher Reynolds number regime (1,700) indicate a dependence on the flow structures

PAGE 117

106 history. In the low Reynolds number regime (100) the viscous effect dissipates the vortex structures quickly and the flow history has almost no impact on the force patterns; symmetric kinematics leads to symmetric force patterns.

PAGE 118

CHAPTER 5 SUMMARY AND FUTURE WORK In this chapter, the progress made regarding the aerodynamics of fixed and flapping wings at low Reynolds number are summarized. Then, future directions of research to obtain a better understanding of unsteady aerodynamics of flexible, fixed and flapping wings are proposed. 5.1 Present Progress of Research Work A computational framework for the fluid-structure interaction problem was developed by coupling a finite volume, pressure-based, Navier-Stokes solver, a finite element structural solver, and a moving grid algorithm to handle shape changes and the resulting grid deformation. The low aspect ratio wing, resulting from the size limitation of MAVs, induces strong tip vortices increase the induced drag and substantially reduce the lift-to-drag ratio. Efforts have been made to numerically probe the effect of endplates in reducing the drag and increasing the lift characteristics of a 15 cm MAV wing. The results demonstrate that the presence of an endplate parallel to the free-stream velocity increases lift without a significant increase in drag. The increased lift is largely realized in the tip region where the endplate reduces the downwash, and increases the effective angle of attack. However, the effectiveness of the endplate diminishes as the angle of attack increases, due to stronger wing tip vortices. The impact of wing flexibility was studied on a 12.7 cm wingspan wing at 6 o angle of attack and Reynolds number of 7.1x10 4 The results showed that the membrane wing 107

PAGE 119

108 configuration vibrates, even in a steady free-stream. Similar experimental (Waszak et al., 2001) and computational (Lian et al., 2003c) investigations on MAVs found that the membrane vibrates at O(100) Hz. Time dependent vertical displacements of the membrane surface are assessed. Although the accuracy of the displacement field is uncertain, as the rigidity of the battens is not precisely accounted for, results obtained from the current approach are consistent with the experimental observations. In terms of aerodynamic forces, for small to moderate angles of attack and before the stall limit, the time-averaged lift and drag coefficients of the membrane wing are comparable to those of a rigid wing. The membrane wing can extend the stall limit over the rigid wing by a substantial margin. The mechanisms of lift generation in flapping flight are investigated for two hovering modes; one is based on the water treading mode and the other on the insect normal hovering mode. The results for the water treading mode with a stroke amplitude of 1 chord and at a Reynolds number of 1,700 are compared with similar numerical simulations. The fluid solver coupled with the moving grid algorithm provides good predictions in terms of force history and average lift coefficient. The numerical results for the normal hovering mode are compared with experimental and prior computational results of Wang et al. (2004). The force history and peak values for a stoke amplitude of 1.4 chords at a Reynolds number of 75 are similar with computational and experimental measurements of Wang et al. (2004). For the normal hovering mode, with a stroke amplitude of 2.4 chords at a Reynolds number of 115, the flow patterns of the present computations share substantial similarities with the experimental flow visualizations of Wang et al. (2004) and further confirm the

PAGE 120

109 capabilities of the present numerical framework to resolve unsteady, low Reynolds number flows. Different lift mechanisms are identified as contributing to the lift generation for two hovering modes at a low Reynolds number of 100. The first mode, based on the water-treading mode, indicates that the sinusoidal-type lift peaks are primarily generated by the rapid pitching-up and delayed-stall mechanism. For the second hovering mode, based on the normal hovering kinematics, the local lift peaks observed at the beginning of the half-strokes indicates (besides rapid pitch-up and delayed-stall mechanisms) that a wake capturing mechanism is present. The average lift coefficient demonstrates that the water treading mode produced more lift than the normal hovering mode for the studied Reynolds number. A variation of the Reynolds number by one order of magnitude significantly changes the force peak values for the hovering mode based on the water treading mode with a large stroke amplitude (1.4 chords). Although the force patterns at both Reynolds numbers share similarity, the stronger convective effect of the higher Reynolds number case exhibit more pronounced asymmetry between forward and backward strokes. 5.2 Future Research Direction Recent numerical and experimental investigations of the flapping flight provided a better understanding of the unsteady lift generation mechanisms. Nevertheless, important work lies ahead to clarify force generation mechanisms and implement the gained knowledge to a working flapping wing MAV. The main research directions proposed to further investigate the aerodynamics of flapping flight are: Three-dimensional flapping wing simulations with comprehensive evaluation of the kinematics (stroke amplitude, rotation modes, phase lag, etc.) and flight parameters (Reynolds number, hovering, forward flight, etc) on the aerodynamic force

PAGE 121

110 generation. It is important to probe the associated flow structures to offer clear insight into the force generation mechanisms. Behavior of the flexible structure of wings undergoing flapping motion and the effect of flexibility on the lift, drag and thrust. Three dimensional interactions of multiple wing configurations (dragonfly) and their efficiency compared with a single wing pair configuration. Parametric evaluation of lift, drag and thrust generation and flapping efficiency for different configurations, leading to optimal selection of kinematics parameters for desired flying characteristics. Modeling of the laminar-to-turbulent transition observed in the Reynolds number range of operation for micro air vehicles.

PAGE 122

LIST OF REFERENCES Albertani R., Hubner P., Ifju P.G., Lind R. and Jackowski J., 2004, Wind Tunnel Testing of Micro Air Vehicles at Low Reynolds Numbers, Paper No. 2004-01-3090, SAE 2004 World Aviation Conference, Reno NV. Alexander D.E., 1986, Wind Tunnel Studies of Turns by Flying Dragonflies, Journal of Experimental Biology, Vol. 122, Issue 1, pp. 81. Anderson Jr. J.D., 1989, Introduction to Flight, 3rd ed. McGraw-Hill, New York. Appa K., 1989, Finite Surface Spline, Journal of Aircraft, Vol. 26, pp. 495, 496. Azuma A., 1992, The Biokinetics of Flying and Swimming, Springer-Verlag, Tokyo. Azuma A., Azuma S., Watanabe I. and Furuta T., 1985, Flight Mechanics of a Dragonfly, Journal of Experimental Biology, Vol. 116, pp: 79-107. Azuma A. and Watanabe T., 1988, Flight Performance of a Dragonfly, Journal of Experimental Biology, Vol. 137, pp. 221. Batina J., 1989, Unsteady Euler Airfoil Solutions Using Unstructured Dynamic Meshes, AIAA Paper 1989-0115. Bendiksen O.O., 1991, A New Approach to Computational Aeroelasticity, AIAA Paper 91-0939-CP. Bennett L., 1970, Insect Flight: Lift and Rate of Change of Incidence, Science, Vol. 167, pp. 177-179. Bennett L., 1977, Clap and Fling Aerodynamics: An Experimental Evaluation, Journal of Experimental Biology, Vol. 69, pp. 261. Betz A., 1912, Ein Beitrag zur Erklrung des Segelfluges, Zeitschrift fr Flugtechnik und Motorluftschiffahrt, Vol. 3, January, pp. 269-272. Birch J.M. and Dickinson M.H., 2001, Spanwise Flow and the Attachment of the Leading-edge Vortex, Nature, Vol. 412, No. 6848, pp. 729-733. Birch J.M. and Dickinson M.H., 2003, The Influence of Wing-wake Interaction on the Production of the Aerodynamic Forces in Flapping Flight, Journal of Experimental Biology, Vol. 206, pp. 2257-2272. 111

PAGE 123

112 Birnbaum W., 1924, Der Schlagflgelpropeller und die kleinen Schwingungen elastisch befestigter Tragfgel, Zeitschrift fr Flugtechnik und Motorluftschiffahrt, Vol. 15, pp. 128-134. Blasius H., 1908, Grenzschichten in Flssigkeiten mit kleiner Reibung, Zeitschrift fr Mathematik und Physik, Vol. 56, pp. 1-37. English translation in NACA TM 1526. Brodsky A.K., 1991, Vortex Formation in the Tethered Flight of the Peacock Butterfly Inachis io L. (Lepidoptera, Nymphalidae) and Some Aspects of Insect Flight Evolution, Journal of Experimental Biology, Vol. 161, November, pp. 77. Buckholz R.H., 1981, Measurements of Unsteady Periodic Forces Generated by the Blowfly Flying in a Wind Tunnel, Journal of Experimental Biology, Vol. 90, pp. 163-173. Byrne D.N., Buchmann S.L. and Spangler H.G., 1988, Relationship Between Wing Loading, Wingbeat Frequency, and Body Mass in Homopterous Insects, Journal of Experimental Biology, Vol. 135, pp. 9-23. Chorin A.J., 1968, Numerical Solution of the Navier-Stokes Equations, Mathematics of Computation, Vol. 22, October, pp. 745-762. Cloupeau, M., Devilliers, J.F. and Devezeaux, D., (1979) Direct Measurements of Instantaneous Lift in Desert Locust; Comparison with Jensens Experiments on Detached Wings, Journal of Experimental Biology, Vol. 80, pp. 1-15. Cooter, R.J. and Baker, P.S., (1977) Weis-Fogh Clap and Fling Mechanism in Locusta, Nature, Vol. 269, pp. 53. Dickinson M.H. and Gtz K.G., 1993, Unsteady Aerodynamic Performance of Model Wings at Low Reynolds Numbers, Journal of Experimental Biology, Vol. 174, January, pp. 45-64. Dickinson M.H. and Gtz K.G., 1996, The Wake Dynamics and Flight Forces of the Fruit Fly, Drosophila melanogaster, Journal of Experimental Biology, Vol. 199, September, pp. 2085-2104. Dickinson M.H., Lehmann F-O. and Gtz K.G., 1993, The Active Control of Wing Rotation by Drosophila, Journal of Experimental Biology, Vol. 182, September, pp. 173-189. Dickinson M.H., Lehmann F-O. and Sane S., 1999, Wing Rotation and the Aerodynamic Basis of Insect Flight, Science, Vol. 284, June, pp. 1954-1960. Ding H., Yang B., Lou M. and Fang H., 2003, New Numerical Method for Two-dimensional Partially Wrinkled Membranes, AIAA Journal, Vol. 41, No. 1, pp. 125.

PAGE 124

113 Ellington C.P., 1978, The Aerodynamics of Normal Hovering Flight: Three Approaches, in Comparative Physiology: Water, Ions and Fluid Mechanics, ed. K. SchmidtNielsen, L. Bolis and S. Maddrell, Cambridge University Press, Cambridge, pp. 327-345. Ellington C.P., 1980, Vortices and Hovering Flight, in Instationare Effekte an Schwingenden Tierflgeln, ed. W. Nachtigall, Franz Steiner, Wiesbaden, pp. 64-101. Ellington C.P., 1984a, The Aerodynamics of Hovering Insect Flight. I. The Quasi-steady Analysis, Philosophical Transactions of the Royal Society of London, Series B, Vol. 305, February, pp. 1-15. Ellington C.P., 1984b, The Aerodynamics of Hovering Insect Flight. II. Morphological Parameters, Philosophical Transactions of the Royal Society of London, Series B, Vol. 305, February, pp. 17-40. Ellington C.P., 1984c, The Aerodynamics of Hovering Insect Flight. III. Kinematics, Philosophical Transactions of the Royal Society of London, Series B, Vol. 305, February, pp. 41-78. Ellington C.P., 1984d, The Aerodynamics of Hovering Insect Flight. V. A Vortex Theory, Philosophical Transactions of the Royal Society of London, Series B, Vol. 305, February, pp. 115-144. Ellington C.P., 1984e, The Aerodynamics of Hovering Insect Flight. IV. Aerodynamic Mechanisms, Philosophical Transactions of the Royal Society of London, Series B, Vol. 305, February, pp. 79-113. Ellington C.P., 1984f, The Aerodynamics of Hovering Insect Flight. VI. Lift and Power Requirements, Philosophical Transactions of the Royal Society of London, Series B, Vol. 305, February, pp. 145-181. Ellington C.P., Berg C. van den, Willmott A.P., Thomas A.L.R., 1996, Leading-edge Vortices in Insect Flight, Nature, Vol. 384, December, pp. 626-630. Eriksson L.E., 1982, Generation of Boundary-conforming Grids Around Wing-body Configurations Using Transfinite Interpolation, AIAA Journal, Vol. 20, No. 10, pp. 1313. Freymuth P., 1990, Thrust Generation by an Airfoil in Hover Mode, Experiments in Fluids, Vol. 9, January, pp. 17-24. Gad-el-Hak M., 2001, Micro-air-vehicles: Can They Be Controlled Better?, Journal of Aircraft, Vol. 38, No. 3, pp. 419.

PAGE 125

114 Gilmanov A. and Sotiropoulos F., 2005, A Hybrid Cartesian/immersed Boundary Method for Simulating Flows with 3D, Geometrically Complex, Moving Bodies, Journal of Computational Physics (Article in Press), January. Goldspink G., 1977, Energy Cost of Locomotion, in: Alexander, R.M., Chapman, G.C. editors, Mechanics and Energetics of Animal Locomotion, Chapman and Hall, London. Gordnier R.E. and Fithen R., 2003, Coupling of a Nonlinear Finite Element Structural Method with a NavierStokes Solver, Computers and Structures, Vol. 81, pp.7589. Gtz K.G., 1987, Course-control, Metabolism and Wing Interference During Ultralong Tethered Flight in Drosophila melanogaster, Journal of Experimental Biology, Vol. 128, March, pp. 35-46. Grasmeyer J.M. and Keennon M.T., 2001, Development of the Black Widow Micro Air Vehicle, AIAA Paper 2001-0127. Greenewalt C.H., 1962, Dimensional Relationships for Flying Animals, Smithsonian miscellaneous collections, Smithsonian Institution, Washington, Vol. 144, No. 2, pp. 1-46. Green A.E. and Adkins J.E., 1960, Large Elastic Deformations. The Clarendon Press, Oxford. Hartwich P.M. and Agrawal S., 1997, Method for Perturbing Multiblock Patched Grids in Aeroelastic and Design Optimization Applications, AIAA Paper 1997-2038. Hill A.V., 1950, The Dimensions of Animals and Their Muscular Dynamics, Science Progress, Vol. 38, pp. 209-230. Ho S., Nassef H., Pornsin-sirirak T.N., Tai Y-C. and Ho C-M., 2003, Unsteady Aerodynamics and Flow Control for Flapping Wing Flyers, Progress in Aerospace Sciences, Vol. 39, pp. 635-681. Ifju P.G., Ettinger S., Jenkins D. and Martinez L., 2001, Composite Materials for Micro Air Vehicles, SAMPE Journal, Vol. 37, pp. 7. Ifju P.G., Jenkins D., Ettinger S., Lian Y., Shyy W. and Waszak R.M., 2002, Flexible-wing-based Micro Air Vehicles, AIAA Paper 2002-0705. Issa R.I., 1985, Solution of the Implicitly Discretised Fluid Flow Equations by Operator-splitting, Journal of Computational Physics, Vol. 62, pp. 40-65. Jackson P.S. and Christie G.W., 1987, Numerical Analysis of Three-dimensional Elastic Membrane Wings, AIAA Journal, Vol. 25, No. 5, pp. 676.

PAGE 126

115 Jameson A., 1991, Time Dependent Calculations Using Multigrid, with Applications to Unsteady Flows Past Airfoils and Wings, AIAA Paper 1991-1596. Jenkins C.H., 1996, Nonlinear Dynamic Response of Membranes: State of the Artupdate, Applied Mechanical Review, Vol. 49, No. 10, pp. S41S48. Jenkins C.H. and Leonard J.W., 1991, Nonlinear Dynamic Response of Membranes: State of the Art, Applied Mechanical Review, Vol. 44, pp.319. Jensen M., 1956, Biology and Physics of Locust Flight. III. The Aerodynamics of Locust Flight, Philosophical Transactions of the Royal Society of London, Series B, Vol. 230, July, pp. 511-552. Jones K.D., Bradshaw C.J., Papadopoulos J. and Platzer M., 2004, Improved Performance and Control of Flapping-wing Propelled Micro Air Vehicles, AIAA Paper 2004-399. Jones K.D. and Platzer M.F., 2003, Experimental Investigation of the Aerodynamic Characteristics of Flapping-wing Micro Air Vehicles, AIAA Paper 2003-0418. Karman T. von and Burgers J.M., 1935, Aerodynamic Theory, Vol. 2, Springer, Berlin. Katzmayr R., 1922, Effect of Periodic Changes of Angle of Attack on Behavior of Airfoils, NACA TM No. 147, October. Kesel A.B., 2000, Aerodynamic Characteristics of Dragonfly Wing Sections Compared with Technical Airfoils, Journal of Experimental Biology, Vol. 203, pp. 3125-3135. Knoller R., 1909, Die Gesetze des Luftwiderstandes, Flugund Motortechnik (Wien), Vol. 3, No. 21, pp. 1-7. LaRoche U. and Palffy S., 1996, Wing Grid, a Novel Device for Reduction of Induced Drag on Wings, Presented at International Council of Aeronautical Sciences (ICAS), Sorrento, Italy. Lehmann F-O., 2000, Flattern fr Flugkrfte, Naturwiss Rundschau, Vol. 623, pp. 223-230. Lehmann F-O., 2004, The Mechanisms of Lift Enhancement in Insect Flight, Naturwissenschaften, Vol. 91, March, pp. 101-122. Lian Y. and Shyy W., 2005, Numerical Simulations of Membrane Wing Aerodynamics for Micro Air Vehicles Applications, Journal of Aircraft, Vol. 42, No. 4, pp. 865-873. Lian Y., 2003a, Membrane and Adaptively-shaped Wings for Micro Air Vehicles, Ph.D. Dissertation, University of Florida, Gainesville, FL.

PAGE 127

116 Lian Y., Shyy W. and Haftka R., 2003b, Shape Optimization of a Membrane Wing for Micro Air Vehicles, AIAA Paper 2003-0106. Lian Y., Shyy W., Viieru D. and Zhang B., 2003c, Membrane Wing Aerodynamics for Micro Air Vehicles, Progress in Aerospace Sciences, Vol. 39, pp. 425-465. Lian Y., Steen J., Trygg-Wilander M. and Shyy W., 2003d, Low Reynolds Number Turbulent Flows Around a Dynamically Shaped Airfoil, Computers and Fluids, Vol. 32, Issue 3, pp. 287. Lian Y., Shyy W., Ifju P.G. and Verron E., 2002, A Computational Model for Coupled Membrane-fluid Dynamics, AIAA Paper 2002-2972. Lighthill M.J., 1973, On the Weis-Fogh Mechanism of Lift Generation, Journal of Fluid Mechanics, Vol. 60, pp. 1. Lilienthal O., 1889, Der Vogelflug als Grundlage der Fliegekunst, ed. R. Gaertners Verlag. Transl. by A.W. Isenthal, 1911, Birdflight as the Basis of Aviation, Longmans, Green & Co., London, UK. Lippisch A.M., 1960, Man-powered Flight in 1929, Journal of Royal Aeronautical Society, Vol. 64, July, pp. 395-398. Lissaman P.B.S., 1983, Low-Reynolds-number Airfoils, Annual Review of Fluid Mechanics, Vol. 15, No. 2, pp. 23. Liu H., Ellington C.P., Kawachi K., Berg C., van den and Willmott A.P., 1998, A Computational Fluid Dynamic Study of Hawkmoth Hovering, Journal of Experimental Biology, Vol. 201, February, pp. 461477. Liu H. and Kawachi K., 1998, A Numerical Study of Insect Flight, Journal of Computational Physics, Vol. 146, October, pp. 124-156. Liu H., Wang H. and Inada Y., 2005, Integrated Modeling and Simulation of Free Flight in Realistic Insect, Presented at the Society of Experimental Biology (SEB) Annual Meeting, July, Barcelona, Spain. Marey E.J., 1868, Determination Experimentale du Mouvement des Ailes des Insectes Pendant le Vol, C. R. Acad. Sci. Paris Vol. 67, pp. 1341. Maxworthy T., 1979, Experiments on the Weis-Fogh Mechanism of Lift Generation by Insects in Hovering Flight. Part 1. Dynamics of the Fling, Journal of Fluid Mechanics, Vol. 93, pp. 47-63. Messiter, A.F., 1970, Boundary Layer Flow Near the Trailing Edge of a Flat Plate, SIAM Journal of Applied Mathematics, Vol. 18, pp. 241-257.

PAGE 128

117 Mittal R., Utturkar Y. and Udaykumar H.S., 2002, Computational Modeling and Analysis of Biomimetic Flight Mechanisms, AIAA Paper 2002-865. Mittal R., Dong H., Bozkurttas M., Von Loebbecke A. and Najjar F., 2006, Analysis of Flying and Swimming in Nature Using an Immersed Boundary Method, AIAA Paper 2006-2867. Mooney, M., 1940, A Theory of Large Elastic Deformation, Journal of Applied Physics, Vol. 11, pp. 582. Mueller T.J., ed. 1989, Proceedings of the Conference on Low Reynolds Number Aerodynamics. University of Notre Dame, Notre Dame, IN. Mueller T.J., ed. 2000, Proceedings of the Conference on Fixed, Flapping and Rotary Wing Vehicles at Very Low Reynolds Numbers. University of Notre Dame, Notre Dame, IN. Mueller T.J. and DeLaurier J.D., 2003, Aerodynamics of Small Vehicles, Annual Review of Fluid Mechanics, Vol. 35, pp. 89. Nachtigall W., 1976, Wing Movements and the Generation of Aerodynamic Forces by Some Medium-sized Insects, in: Rainey, R.C., (ed) Insect flight. Blackwell, Oxford, pp. 31-49. Nelsen J.N., 1963, Theory of Flexible Aerodynamics Surfaces, Journal of Applied Mechanics Vol. 30, No. 4, pp. 35. Newman B.G., 1987, Aerodynamics Theory for Membrane and Sails, Progress in Aerospace Sciences, Vol. 24, pp.1. Norberg R.A., 1975, Hovering Flight of the Dragonfly Aeschna juncea L, in: Wu, T.Y., Brokaw, C.J., Brennen, C., (eds) Kinematics and Aerodynamics, Vol 2, Plenum, New York, pp 763. Oden J.T. and Sato T., 1967, Finite Strains and Displacements of Elastic Membranes by the Finite Element Method, International Journal of Solids and Structures, Vol. 3, pp. 471. Oliveira P.J. and Issa R.I., 2001, An Implicit PISO Algorithm for the Computation of Buoyancy-driven Flows, Numerical Heat Transfer Part B, Vol. 40, Issue 6, pp.473. Osborne, M.F.M., 1951, Aerodynamics of Flapping Flight with Applications to Insects, Journal of Experimental Biology, Vol. 28, pp. 221-245. Patankar S.V. and Spalding D.B., 1972, A Calculation Procedure for Heat, Mass and Momentum Transfer in Three-dimensional Parabolic Flows, International Journal of Heat and Mass Transfer Vol.15, pp.1787.

PAGE 129

118 Pennycuick C.J., 1996, Wingbeat Frequency of Birds in Steady Cruising Flight: New Data and Improved Predictions, Journal of Experimental Biology, Vol. 199, pp.1613-1618. Pornsin-Sirirak T.N., Lee S.W., Nassef H., Grasmeyer J., Tai Y.C., Ho C.M. and Keennon M., 2000, MEMS Wing Technology for a Battery-powered Ornithopter, Micro Electro Mechanical Systems, 2000, The 13 th Annual International Conference on, Miyazaki, Japan, 23-27 January, pp. 799-804. Prandtl L., 1904, ber Flssigkeitsbewegung bei sehr kleiner Reibung, Proceedings of the Third International Mathematics Congress, Heidelberg 1904. English translation in NACA TM 452 (1928). Prandtl L., 1924, Uber die Entstehung von Wirbeln in einer idealen Flussigkeit, Vortrage zur Hydround Aerodynamik, Berlin, p. 5. Pulliam T., 1993, Time Accuracy and the Use of Implicit Methods, AIAA Paper 93-3360-CP. Ramamurti R. and Sandberg W., 1994, Evaluation of a Scalable 3-D Finite Element Incompressible Flow Solver, AIAA Paper 1994-0756. Ramamurti R. and Sandberg W., 2002, A Three-dimensional Computational Study of the Aerodynamic Mechanisms of Insect Flight, Journal of Experimental Biology, Vol. 205, May, pp. 1507-1518. Rayner J.M.V., 1979a, A Vortex Theory of Animal Flight. Part 1. The Vortex Wake of a Hovering Animal, Journal of Fluid Mechanics, Vol. 91, pp. 697-730. Rayner J.M.V., 1979b, A Vortex Theory of Animal Flight. Part 2. The Forward Flight of Birds, Journal of Fluid Mechanics, Vol. 91, pp. 731-763. Reavis M.A., and Luttges M.W., 1998, Aerodynamic Forces Produced by a Dragonfly, AIAA Paper 1988-330. Robinson B.A., Yang H.T.Y. and Batina J.T., 1991, Aeroelastic Analysis of Wings Using the Euler Equations with a Deforming Mesh, Journal of Aircraft, Vol. 28, No. 11, pp. 781. Rogers S.E. and Kwak D., 1990, Upwind Differencing Scheme for the Time-accurate Incompressible NavierStokes Equations, AIAA Journal Vol. 28, February, pp. 253-262. Roger S.E., Kwak D. and Kiris C., 1991, Steady and Unsteady Solutions of the Incompressible Navier-Stokes Equations, AIAA Journal, Vol. 29, pp. 603-610. Rozhdestvensky K.V. and Ryzhov V.A., 2003, Aerohydrodynamics of Flapping-wing Propulsors, Progress in Aerospace Sciences, Vol. 39, November, pp. 585-681.

PAGE 130

119 Rumsey C.L., Sanetrik M.D., Biedron R.T., Melson N.D. and Parlette E.B., 1996, Efficiency and Accuracy of Time-accurate Turbulent NavierStokes Computations, Computers and Fluids, Vol. 25, Issue 2, pp. 217-236. Sane S.P., 2003, The Aerodynamics of Insect Flight, Journal of Experimental Biology, Vol. 206, pp. 4191-4208. Sane S.P. and Dickison M.H., 2002, The Aerodynamic Effects of Wing Rotation and a Revised Quasi-steady Model of Flapping Flight, Journal of Experimental Biology, Vol. 205, pp. 1087-1096. Saharon D. and Luttges M. W., 1987, Three-dimensional Flow Produced by a Pitching-plunging Model Dragonfly Wing, AIAA Paper 1987-0121. Saharon D. and Luttges M.W., 1988, Visualization of Unsteady Separated Flow Produced by Mechanically Driven Dragonfly Wing Kinematics Model, AIAA Paper 1988-0569. Savage S.B., Newman B.G. and Wong D.T.M., 1979, The Role of Vortices and Unsteady Effects During the Hovering Flight of Dragonflies, Journal of Experimental Biology, Vol. 83, pp. 59-77. Schlichting H., 1979, Boundary-layer Theory, 7th ed. McGraw-Hill, New York. Schmidt-Nielsen K., 1972, Locomotion: Energy Cost of Swimming, Flying and Running, Science, Vol. 177, pp. 222-228. Schuster D., Vadyak J. and Atta E., 1990, Static Aeroelastic Analysis of Fighter Aircraft Using a Three-dimensional NavierStokes Algorithm, Journal of Aircraft, Vol. 27 No. 9, pp. 820. Shyy W., 1994, Computational Modeling for Fluid Flow and Interfacial Transport. Elsevier, Amsterdam, The Netherlands. Shyy W., Berg M. and Ljungqvist D, 1999, Flapping and Flexible Wings for Biological and Micro Air Vehicles, Progress in Aerospace Sciences, Vol. 35, July, pp. 455506. Shyy W. and Mittal R., 1998, Solution Methods for the Incompressible NavierStokes Equations, In: Johnson R, ed. Handbook of Fluid Dynamics. CRC Press, Boca Raton, FL, pp. 31-1 to 31-33. Shyy W., Udaykumar H.S., Rao M.M. and Smith R.W., 1996, Computational Fluid Dynamics with Moving Boundaries, Taylor & Francis, Washington, DC. Smith M.J.C., 1996, Simulating Moth Wing Aerodynamics Towards the Development of Flapping-wing Technology, AIAA Journal, Vol. 34, No. 7, pp. 1348-1355.

PAGE 131

120 Smith M.J., Hodges D.H. and Cesnik C.E.S., 2000, Evaluation of Computational Algorithms Suitable for Fluidstructure Interactions, Journal of Aircraft, Vol. 37, No. 2, pp. 282. Smith R.W. and Shyy W., 1997, Incremental Potential Flow Based Membrane Wing Element, AIAA Journal, Vol. 35, No 5, pp. 782. Spedding G.R. and Maxworthy T., 1986, The Generation of Circulation and Lift in a Rigid Two-dimensional Fling, Journal of Fluid Mechanics, Vol. 165, pp. 247-272. Srygley R.B. and Thomas A.L.R., 2002, Unconventional Lift-generating Mechanisms in Free-flying Butterflies, Nature, Vol. 420, No. 6916, pp. 660-664. Stewartson K., 1974, Multistructured Boundary Layers on Flat Plates and Related Bodies, Advanced Applied Mechanics, Vol. 14, Academic Press, New York, pp. 146-239. Stokes G.G., 1851, On the Effect of the Internal Friction of Fluids on the Motion of Pendulums, Cambridge Philosophical Transactions, Vol. 9, pp. 8-106. Sun M. and Tang J., 2002a, Unsteady Aerodynamic Force Generation by a Model Fruit Fly Wing in Flapping Motion, Journal of Experimental Biology, Vol. 205, January, pp. 55-70. Sun M. and Tang J., 2002b, Lift and Power Requirements of Hovering Flight in Drosophila virilis, Journal of Experimental Biology, Vol. 205, August, pp. 2413-2427. Sun M. and Lan S., 2004, A Computational Study of the Aerodynamic Forces and Power Requirements of Dragonfly (Aeschna juncea) Hovering, Journal of Experimental Biology, Vol. 207, May, pp. 1887. Sunada S., Kawachi K., Watanabe I. and Azuma A., 1993, Fundamental Analysis of Three-dimensional Near Fling, Journal of Experimental Biology, Vol. 183, October, pp. 217. Tang J. and Sun M., 2001, Force and Flow Structures of a Wing Performing Flapping Motion at Low Reynolds Number, Acta Mechanica, Vol. 152, No. 1-4, pp. 35-48. Tang J., Viieru D. and Shyy W., 2006, Effects of Reynolds Number and Stroke Amplitude on Hovering Aerodynamics, submitted for review for the 45th AIAA Aerospace Sciences Meeting and Exhibit, Reno NV, Jan. 2007. Tani I., 1964, Low-speed Flows Involving Bubble Separations, In: Kuchemann, D., Sterne, H.G., editors. Progress in Aeronautical Science, Vol. 5, Pergamon, New York, pp. 70.

PAGE 132

121 Tannehill J.C., Anderson D.A. and Pletcher R.H., 1997, Computational Fluid Mechanics and Heat Transfer, 2nd ed., Taylor and Francis, Bristol, PA. Thakur S., Wright J. and Shyy W., 2002, STREAM: A Computational Fluid Dynamics and Heat Transfer NavierStokes Solver. Theory and Applications, Streamline Numerics, Inc., and Computational Thermo-Fluids Laboratory, Department of Mechanical and Aerospace Engineering Technical Report, Gainesville, FL. Thomas P.D. and Lombard, C.K., 1979, Geometric Conservation Law and Its Application to Flow Computation on Moving Grids, AIAA Journal Vol. 17, No. 10, pp. 1030-1037. Thwaites B., 1961, Aerodynamics Theory of SailsPart I, Proceedings of the Royal Society of London, Series A, Vol. 261, pp.402. Torres G.E. and Mueller T.J., 2001, Aerodynamic Characteristics of Low Aspect Ratio Wings at Low Reynolds Number, In: Mueller, T.J., ed. Fixed and Flapping Wing Aerodynamics for Micro Air Vehicle Applications, vol. 195, AIAA, Reston, VA, pp. 191. Udaykumar H.S., Mittal R. and Shyy W., 1999, Computation of Solid-liquid Phase Fronts in the Sharp Interface Limit on Fixed Grids, Journal of Computational Physics, Vol. 153, August, pp. 535-574. Udaykumar H.S., Mittal R., Rampunggoon P. and Khanna A., 2001, A Sharp Interface Cartesian Grid Method for Simulating Flows with Complex Moving Boundaries, Journal of Computational Physics, Vol. 174, November, pp. 345. Van Doormaal J.P. and Raithby G.D., 1985, Enhancements of the SIMPLE Method for Predicting Incompressible Fluid Flows, Numerical Heat Transfer Vol.67, pp.147163. Verron E., Marckmann G. and Peseux B., 2001, Dynamic Inflation of Non-linear Elastic and Viscoelastic Rubber-like Membranes, International Journal for Numerical Methods in Engineering, Vol. 50, Issue 5, pp. 1233. Viieru D., Lian Y., Shyy W. and Ifju P.G., 2003, Investigation of Tip Vortex on Aerodynamic Performance of a Micro Aerial Vehicle, AIAA Paper 2003-3597. Viieru D., Albertani R., Shyy W. and Ifju P.G., 2005, Effect of Tip Vortex on Wing Aerodynamics of Micro Aerial Vehicles, Journal of Aircraft, Vol. 42, No. 6, pp. 1530-1536. Wagner H., 1925, ber die Entstehung des dynamischen Auftriebes von Tragflgeln, Z. Angew. Math. Mech. Vol. 5, pp. 17-35. Walker A.J., 2002, Rotational Lift: Something Different or More of the Same? Journal of Experimental Biology, Vol. 205, pp. 3783.

PAGE 133

122 Walker G.T., 1925, The Flapping Flight of Birds. I., Journal of the Royal Aeronautical Society Vol. 29, pp. 590. Walker G.T., 1927, The Flapping Flight of Birds. II., Journal of the Royal Aeronautical Society Vol. 31, pp. 337. Walker P.B., 1931, Experiments on the Growth of Circulation About a Wing and an Apparatus for Measuring Fluid Motion, Aeronautical Research Committee (Great Britain) Technical Report No. 1402. Wang Z.J., 2000a, Two Dimensional Mechanism for Insect Hovering, Physical Review Letters, Vol. 85, September, pp. 2216-2219. Wang Z.J., 2000b, Vortex Shedding and Frequency Selection in Flapping Flight, Journal of Fluid Mechanics, Vol. 410, May, pp. 323-341. Wang Z.J., 2005, Dissecting Insect Flight, Annual Review of Fluid Mechanics, Vol. 37, January, pp. 183-210. Wang Z.J., Birch J.M. and Dickinson M.H., 2004, Unsteady Forces and Flows in Low Reynolds Number Hovering Flight: Two-dimensional Computations vs Robotic Wing Experiments, Journal of Experimental Biology, Vol. 207, January, pp. 449-460. Waszak R. M., Jenkins N.L. and Ifju P.G., 2001, Stability and Control Properties of an Aeroelastic Fixed Wing Micro Aerial Vehicle, AIAA Paper 2001-4005. Weis-Fogh T., 1973, Quick Estimate of Flight Fitness in Hovering Animals, Including Novel Mechanisms for Lift Production, Journal of Experimental Biology, Vol. 59, pp. 169-230. White, F.M., 1991, Viscous Fluid Flow, 2nd ed. McGraw-Hill, New York. Wilkin P.J. and Williams M.H., 1993, Comparison of the Instantaneous Aerodynamic Forces on a Sphingid Moth with Those Predicted by Quasi-steady Aerodynamic Theory, Physiol. Zool., Vol. 66, pp. 1015-1044. Willmott A.P. and Ellington C.P., 1997a, The Mechanics of Flight in the Hawkmoth Manduca Sexta. I. Kinematics of Hovering and Forward Flight, Journal of Experimental Biology, Vol. 200, November, pp. 2705-2722. Willmott A.P. and Ellington, C.P., 1997b, Measuring the Angle of Attack of Beating Insect Wings: Robust Three-dimensional Reconstruction from Two-dimensional Images, Journal of Experimental Biology, Vol. 200, pp. 2693-2704. Wilson E.L., Farhoomand I. and Bathe K.J., 1973, Nonlinear Dynamic Analysis of Complex Structures, Earthquake Engineering and Structures Dynamics, Vol. 1, pp. 241.

PAGE 134

123 Ye T., Mittal R., Udaykumar H. S. and Shyy W., 1999, An Accurate Cartesian Grid Method for Simulation of Viscous Incompressible Flows with Complex Immersed Boundaries, Journal of Compuational Physics, Vol. 156, December, pp. 209. Zbikowski R., 2002, On Aerodynamic Modeling of an Insect-like Flapping Wing in Hover for Micro Air Vehicles, Philosophical Transactions of the Royal Society of London, Series A, Vol. 360, February, pp. 273-290. Zwaan R.J. and Prananta B.B., 2002, Fluid/structure Interaction in Numerical Aeroelastic Simulation, International Journal of Non-Linear Mechanics, Vol. 37, pp. 987.

PAGE 135

BIOGRAPHICAL SKETCH Dragos Viieru was born in Ploiesti, Romania, in December 1975. He received his B.S. degree in aerospace engineering from the Politehnica University of Bucharest, Romania, in July 1999. In June 2000 he got his Advanced Studies Diploma from the same department. In October 1999 he joined the Romanian Civil Aeronautic Authority as a Flight Operations Inspector and in 2001 he was appointed as Aerodynamic Examiner. In spring 2002 he started his Ph.D. studies in the Department of Mechanical and Aerospace Engineering, University of Florida. His research interests include fluid-structure interaction, aeroelasticity and computational fluid dynamics (CFD). 124


Permanent Link: http://ufdc.ufl.edu/UFE0015704/00001

Material Information

Title: Flapping and Fixed Wing Aerodynamics of Low Reynolds Number Flight Vehicles
Physical Description: Mixed Material
Copyright Date: 2008

Record Information

Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
System ID: UFE0015704:00001

Permanent Link: http://ufdc.ufl.edu/UFE0015704/00001

Material Information

Title: Flapping and Fixed Wing Aerodynamics of Low Reynolds Number Flight Vehicles
Physical Description: Mixed Material
Copyright Date: 2008

Record Information

Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
System ID: UFE0015704:00001


This item has the following downloads:


Full Text












FLAPPING AND FIXED WING AERODYNAMICS
OF LOW REYNOLDS NUMBER FLIGHT VEHICLES













By

DRAGOS VIIERU


A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY

UNIVERSITY OF FLORIDA


2006

































Copyright 2004

by

Dragos Viieru















ACKNOWLEDGMENTS

I would like to express my deepest appreciation to my advisor, Dr. Wei Shyy, for

his guidance and support. I am also grateful for his patience and constant encouragement

throughout my research activity. I owe special thanks to Dr. Peter Ifju and his micro air

vehicles research group for sharing their knowledge and experience. I would also like to

thank Dr. Renwei Mei, Dr. Corin Segal and Dr. William Hager for serving as members of

my supervisory committee.

I acknowledge the close collaboration with Dr. Roberto Albertani regarding

experimental aspects of micro air vehicle aerodynamics and Kyuho Lee who kindly

provided me with his wing shape designs that made part of the work in chapter 3

possible. I would like to thank all my colleagues of the computational thermo-fluids

group for their friendship and support, in particular, Dr. Siddharth Thakur and Dr. Ramji

Kamakoti for generously sharing their time to answer my queries. I would like to thank

Dr. Yongsheng Lian for his collaboration during my first year at University of Florida.

My thanks also go to Dr. Jian Tang, who shared his broad experience on flapping flight

aerodynamics.

Finally, I would like to thank to my parents for their understanding, patience and

encouragement when it was most needed.
















TABLE OF CONTENTS



A C K N O W L E D G M E N T S ................................................................................................. iii

LIST OF TABLES ........... .............................. ................. vi

L IST O F F IG U R E S .... ...... ................................................ .. .. ..... .............. vii

A B ST R A C T .......... ..... ...................................................................................... x


CHAPTER

1. INTRODUCTION TO FLAPPING FLIGHT ................................... ...............1

1 .1 In tro d u ctio n ............................................................................... 1
1.2 F lapping W ing in N ature ............................................................. .....................2
1.2.1 W ing M options ....................................... ........... .... .. ........ ..
1.2 .2 F light M odes in B irds ............................................................................ .
1.2.3 Beating Frequency and Scaling ..................................... ...............5
1.3 Insect Flight ............................10..............................................10
1.3.1 General Characteristics...................... ............. ....... ...............10
1.3.2 Experimental Efforts in Studying the Insect Flight Mechanisms .............11
1.3.3 Computational Efforts in Studying the Insect Flight Mechanisms ............13
1.4 Unsteady mechanisms identified to enhance lift in insects .................................18
1.4.1 Analytical and quasi-steady models of insect flight...............................18
1.4 .2 W agner E effect .................................................................. .............. ..20
1.4.3 Leading-edge Vortex and Delayed Stall................. ............................21
1.4.4 Rotational Forces ...................................... .............. ................24
1.4.5 W ing-w ake Interaction ................................................... .................. 25
1.4.6 Clap and Fling M echanism ..................................... ......... ............... 27
1.5 Objective and O utline ...................................................................... 28

2. COMPUTATIONAL METHODS FOR COUPLED FLUID-STRUCTURE
IN TERA C TION S ............................................ .. .. ........... ......... 30

2.1 O overview ....................................................... ..... ........ ..............30
2.2 Governing Equations for Fluid Flow s ...................................... ............... 31
2 .3 F lu id F low S olv ers........... ..... ......................................................... .... .... .... .. 33
2.4 Geometric Conservation Law ...................................... ......... ............... 34









2.5 M oving G rid Technique ............................................... ............................ 35
2.5.1 M aster-slave C concept ........................................................... .................... 35
2.5.2 Grid Re-arrangement Algorithm with Weighting Functions..................39
2.5.3 Grid Re-arrangement Algorithm Based on Spring Analogy Concept........41
2.6 M em brane Structural Solver........................................... .......................... 42
2.6.1 M ooney-Rivilin M odel ..................................................................44
2.6.2 Solution of the Dynamic Equations................................... .............45
2.7 Coupling between Fluid Solver and Structural Solver .............. ...................47

3. FIXED WING MICRO AIR VEHICLES ...................................... ............... 51

3 .1 O v erv iew ................... ...................................... .......... ................ 5 1
3.2 Rigid W ing M icro A ir V vehicles ........................................ ........ ............... 54
3.2.1 W ing G eom etry ...................... ................ .................. ........ 54
3.2.2 Computational Setup .............................. ...............54
3.2.3 Tip Vortex Effect on the Rigid Wing Aerodynamics.............................. 56
3.3 Flexible W ing M icro Air Vehicles ................. ........................... ............... 65
3.4 Summary and Conclusions for Fixed Wing Aerodynamics ..............................70

4. AERODYNAMICS OF LOW REYNOLDS NUMBER HOVERING AIRFOILS...72

4.1 V alidation of Viscous Force Com putation .................................... ....................72
4.1.1 Steady Flow over a Flat Plate at Zero Angle of Attack..............................73
4.1.2 Flow over an Oscillating Flat Plate in a Quiescent Environment ..............76
4.2 Kinematics of 2-D Airfoil in Hovering Flight..............................................79
4.2.1 Kinematics of "W ater Treading" M ode ...................................................80
4.2.2 Kinematics of "Normal" Hovering Mode ..........................................82
4.3 "Water Treading" Hovering Mode ..................... ..... ...............82
4.3.1 Grid and Time Increment Sensitivity Analysis .......................................84
4.3.2 Com prison w ith Sim ilar W ork......................................................... ........ 87
4.3.3 Lift Generation Mechanism for "Water Treading" Hovering Mode..........87
4.4 "N orm al" H covering M ode .............................................................................91
4.4.1 G rid R efinem ent Study.......................................................... ..................91
4.4.2 Comparison with Similar Computational and Experimental Results.........93
4.5 Lift Generation Mechanisms in "Normal" and "Water Treading" Hovering
M o d e s .............................. .... ... .. ......... ......... .. ............. ............. 9 6
4.6 Summary and Conclusion for the Aerodynamics of Hovering Airfoils.............104

5. SUMMARY AND FUTURE WORK ........................................... ............... 107

5.1 Present Progress of Research Work................................................. 107
5.2 Future Research Direction .................................................................... 109

L IST O F R E F E R E N C E S ........................ .. ........... .............................. .....................

BIOGRAPHICAL SKETCH .......................... ....... ..................... ................... 124
















LIST OF TABLES


Table page

1-1. Variety of beating frequency and Reynolds number for birds and insects. ............7

1-2. Experimental robotic and mechanic devices for insect flapping wing study..........14

1-3. Principal numerical studies of insect flapping flight................... ...... .. .........18

2-1. Moving boundaries techniques. Advantages and disadvantages ..........................42

3-1. Aerodynamic performances for flexible and rigid wing. ......................................69

4-1. Numerical and analytical friction coefficient values for flat plate........................75

4-2. Grid and time increment details for the "water treading" hovering cases .............84

4-3. Grid and time increment details for the "normal" hovering cases ........................91

4-4. Kinematics parameters for the experimental and computational setup...................94

4-5. Kinematic parameters for "water treading" and "normal" hovering modes. The
Reynolds number for both cases is 100.............. ............ ............................... 97

4-6. Parameters for "water treading" hovering mode employed for Reynolds number
effect stu dy ..................................................... ................. 10 1
















LIST OF FIGURES


Figure p

1-1. W ing beating m otions.. .............................. ... .........................................

1-2. H covering flight. .......... ... .... .. .. ..... ..... ...... .... ...... .......... .... ....

1-3. Wing beat frequency versus the body mass for different insect species ................9

1-4. Cross-sectional corrugations of dragonfly wings .................................................12

1-5. Leading edge vortex on a hawkmoth wing. .................................. ...............15

1-6. Flow around a thin airfoil...............................................................................22

1-7. Leading edge vortex (LEV) structures at different Reynolds numbers...................24

1-8. Momentum transfer due to wing-wake interaction..............................................26

1-9. Schematic of clap and fling mechanism................... ............. ............... 27

2-1. Illustration of the moving grid technique .............. ........ .................38

2-2. Master-slave algorithm performance for large rotations................................. 39

2-3. Current node neighbors position for grid re-arrangement algorithm ....................40

2-4. Grid quality for the grid re-arrangement algorithm............... ....................... 41

2-5. Generic algorithm for fluid-structure interaction .......................................... 49

3-1. Fixed and flapping wing M AVs. ........................................ ......................... 53

3-2. M A V 's w ing shapes. .............................................................................. .............54

3-3. Computational domain for the MAV wing. ................................. .................55

3-4. Streamlines, vorticity and pressure on the 6" MAV wing's upper surface as a
function of the angle of attack. ...................................................... ..................57

3-5. W ing shape geom etry ...................................................................... .................. 59









3-6. Vorticity magnitude contours behind the wing at 6 angle of attack. ...................60

3-7. Pressure contours on the wing's upper surface and streamlines slightly above the
w ing at 6 angle of attack.. .................................... .......................... .....................61

3-8. Pressure contours on the wing's lower surface and streamlines slightly below the
w ing at 6 angle of attack. .............................................. ............................. 62

3-9. Pressure coefficient on the wing surface at x/c = 0.34 and 6 angle of attack. .......63

3-10. Pressure coefficient on the wing surface at x/c = 0.53 and 6 angle of attack........63

3-11. Span-wise lift distribution at 6 angle of attack. .............................................. 64

3-12. Span-wise drag distribution at 6 angle of attack. ................................................65

3-13. Surface grid for the 5" M AV wing............................................................ ........ 66

3-14. History of the vertical displacements of two points, one near the wing tip and the
other one between the two battens. .............................................. ............... 67

3-15. Vertical displacement contours at various times .................................... .........68

3-16. Aerodynamic forces history for the 5" flexible wing MAV. ................................68

4-1. Numerical and analytical velocity distribution in the boundary layer along different
locations on a flat plate and different Reynolds numbers. ......................................74

4-2. Numerical and analytical friction coefficient Cffor a flat plate versus Reynolds
num ber ......................... ......................... ........................ 76

4-3. Skin-friction and wall velocity for an oscillating plate. ........................................77

4-4. Velocity distribution above an oscillating plate............... ............................... 78

4-5. Hummingbird in hovering flight. ........................................ ........................ 79

4-6. The schem atics of tw o hovering styles .............. ............. ................. ............. 80

4-7. Schematics of the "water treading" hovering mode ............... ..... ............... 81

4-8. Schematics of "normal" hovering mode. .................................. ....... ............ 82

4-9. The coarse multi-block grid used in the current work. The grid has 90 points
around the airfoil and 70 points in the radial direction. ........................................83

4-10. Lift coefficient for three periods and different grid sizes and 6t=0.01....................85

4-11. Drag coefficient for three periods and different grid sizes and 6t=0.01..................85









4-12. Lift coefficient for three periods and different time increments ...........................86

4-13. Drag coefficient for three periods and different time increments ..........................86

4-14. Lift and Drag coefficients for one period ..................................... ............... 88

4-15. Kinematics of "water treading" hovering mode and the aerodynamic forces for one
com p lete strok e................................................. ................ 89

4-16. Flow field snapshots of the "water treading" hovering mode...............................90

4-17. Lift coefficient of the "normal" hovering mode for three periods and different grid
sizes and St=0.01. ......................................................................92

4-18. Lift coefficient history for the "normal" hovering mode. .....................................95

4-19. Vorticity plot for "normal hovering" mode with h/c=2.4 and Re= 115...............96

4-20. One cycle force history for two hovering modes.. ............................................... 98

4-21. Vorticity contours for two hovering modes.. ...................................................99

4-22. Lift coefficient for the "water treading" mode ............................................... 102

4-23. Vorticity contours for the "water treading" mode...................................103

4-24. Pressure distribution on the airfoil surface for the "water treading" mode...........104















Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy

FLAPPING AND FIXED WING AERODYNAMICS
OF LOW REYNOLDS NUMBER FLIGHT VEHICLES

By

Dragos Viieru

August 2006

Chair: Wei Shyy
Major Department: Mechanical and Aerospace Engineering

Lately, micro air vehicles (MAVs), with a maximum dimension of 15 cm and

nominal flight speed around 10m/s, have attracted interest from scientific and engineering

communities due to their potential to perform desirable flight missions and exhibit

unconventional aerodynamics, control, and structural characteristics, compared to larger

flight vehicles. Since MAVs operate at a Reynolds number of 105 or lower, the lift-to-

drag ratio is noticeably lower than the larger manned flight vehicles. The light weight and

low flight speed cause MAVs to be sensitive to wind gusts. The MAV's small overall

dimensions result in low aspect ratio wings with strong wing tip vortices that further

complicate the aerodynamics of such vehicles. In this work, two vehicle concepts are

considered, namely, fixed wings with flexible structure aimed at passive shape control,

and flapping wings aimed at enhancing aerodynamic performance using unsteady flow

fields. A finite volume, pressure-based Navier-Stokes solver along with moving grid

algorithms is employed to simulate the flow field. The coupled fluid-structural dynamics









of the flexible wing is treated using a hyperelastic finite element structural model, the

above-mentioned fluid solver via the moving grid technique, and the geometric

conservation law. Three dimensional aerodynamics around a low aspect ratio wing for

both rigid and flexible structures and fluid-structure interactions for flexible structures

have been investigated. In the Reynolds numbers range of 7x104 to 9x104, the flexible

wing exhibits self-initiated vibrations even in steady free-stream, and is found to have a

similar performance to the identical rigid wing for modest angles of attack. For flapping

wings, efforts are made to improve our understanding of the unsteady fluid physics

related to the lift generation mechanism at low Reynolds numbers (75 to 1,700).

Alternative moving grid algorithms, capable of handling the large movements of the

boundaries (characteristic of flapping wing kinematics) are tested. Two main hovering

modes are investigated and compared with experimental and other computational efforts.

The analysis shows that, while delayed-stall and rapid pitch-up mechanisms are

responsible for most of the lift generation at a Reynolds numbers of 0(100) and stroke

amplitudes of 0(1 chord), other mechanisms, including wake-capturing, are identified to

contribute to the overall lift/drag force generation. The effect of the Reynolds number on

hovering airfoil aerodynamics is also probed.














CHAPTER 1
INTRODUCTION TO FLAPPING FLIGHT

1.1 Introduction

From early times, the idea of utilizing the thrust generated from flapping wings for

man-made flying vehicles emerged from observations of birds. One of the first attempts

to explain and implement the mechanism of thrust generated from flapping wings to a

flying vehicle was made by Leonardo da Vinci in 1490. Numerous attempts to develop

flying vehicles based on flapping wing aerodynamics were made at the end of the 19th

century and the beginning of the 20th century (Lilienthal, 1889; Lippisch, 1960), but

most of them failed due to inadequate scientific knowledge.

One of the first attempts to explain the physics of flapping wings was made by

Knoller (1909) and Betz (1912). They noted that the wing oscillations produce both thrust

and lift components of the aerodynamic force. In 1922, Katzmayr experimentally

measured the thrust generated by a rigid wing in oscillating flow, confirming the Knoller-

Betz effect. Based on Prandtl's formulation (Prandtl, 1924), Birnbaum (1924) developed

one of the first analytical approaches for calculating thrust and efficiencies of flat plates

in plunging and pitching motion at moderate reduced frequencies. Karman and Burgers

(1935) suggested an explanation for drag or thrust generation of oscillating airfoils based

on the position and orientation of the vortices shed in the wake.

During the second half of the 20th century, research on flapping wing propulsion

intensified, fueled by the fact that numerous biological organisms using flapping wing

propulsion (in air or in water) possess superior aero-hydrodynamic properties and use









energy more efficiently compared with those of the man-made flying vehicles, surface

and underwater ships. Even with the rapid advance of aeronautical knowledge in the last

100 years, insects and birds still amaze us with their impressive performances. For

example, a supersonic aircraft such as the SR-71 Blackbird flying nearly 2000 mph

covers about 32 body lengths per second while a common pigeon flying at 50 mph covers

nearly 75 body lengths (Shyy et al., 1999).

Therefore, nature's flapping wing systems are a precious source of inspiration for

efficient propulsion and lift generation devices used in micro air vehicles. The use of

flapping wing devices is very attractive, because they are able to operate in different

regimes of motion, combine the function of control device, stabilizer and propulsor, and

provide high maneuverability. However, the development of micro air vehicles by

directly copying certain properties of flying organisms is unlikely to be successful. It is

important to thoroughly analyze the applicability of the solutions found by nature for

technical use, and this requires a multi-disciplinary effort involving aero-hydrodynamics,

advanced materials, aero/hydro/servo/elasticity, micro-electro-mechanical systems

(MEMS) technologies, etc.

1.2 Flapping Wing in Nature

If we define energy cost as the energy required to move one kilogram a distance of

one kilometer, and plot it against the body mass of birds, terrestrial animals, and fish

(Schmidt-Nielsen, 1972; Goldsping, 1977), one can see that for a given body size,

running is the most expensive in terms of energy costs, while flying and swimming are

the cheapest. This may be the main reason why nearly a million insect species and more

than 9,000 warm blooded animals have taken to the skies. To counter the gravity force









and propel themselves against the aerodynamic drag, these creatures exclusively use the

beating motion of wings.

Different aspects of flapping wing flight for both birds and insects, such as type of

flight, beating frequency, and scaling, are presented next.

1.2.1 Wing Motions

Azuma (1992) identified four fundamental motions of a beating wing (Fig. 1-1):

1. an out-of-plane motion called "flapping" (flapping angle w );

2. an in-plane motion called "lagging" (lag or azimuthal angle );

3. a twisting motion of the wing pitch called "feathering" (feathering angle 0);

4. an alternatively extending and contracting motion of the wingspan called
"spanning."

The stroke plane is defined by the elevation angle, ,. The feathering motion is

also called supinationn" for positive pitch and pronationn" for negative pitch.


Wing element


Horizontal axis


Vertical axis


Figure 1-1. Wing beating motions. Sketched from Azuma (1992).









Insects employ the same fundamental wing motions, except for spanning, since the

structure of the insect wings do not allow it.

1.2.2 Flight Modes in Birds

Birds can adopt different flight modes to fulfill their objectives in flight by

adjusting the configuration of their wings or wing profile. This allows them to obtain the

optimal wing configuration for each flight mode and each phase of wing beat (Azuma,

1992). In general birds use the following types of flight (Azuma, 1992):

Cruising flight. Cruising flight is a steady level flight in which the almost fully

extended wing acts to give lift force at the inner part of the wing and thrust at the outer

part of the wing with minimum power consumption.

Gliding or soaring flight. Gliding or soaring flight is performed without beating the

wings, except while maneuvering. The wings are fully extended in order to increase the

wings' aspect ratio and obtain a minimum gliding angle (max CL/CD) or a minimum rate

of descent (max C /2 ICD). Here CL and CD are the lift and drag coefficient, respectively.

Diving flight. Diving flight is a steep descent used for preying or landing. The

wings are partially retracted to reduce drag. The birds control the speed and direction

with the help of the tail.

Hovering flight. Hovering flight is the flight mode in which birds can temporarily

maintain their position in calm air. During hovering, the bird wings move in a complex

pattern, with the wing tip describing an elliptical orbit. During the downstroke, the wings

are extended to provide more lift, while on the upstroke the wings are flexed backward to

reduce drag. This is usually called "avian stroke" (Azuma, 1992) and is illustrated in Fig.

1-2(A). Hummingbirds and insects can sustain longer periods of hovering flight by









beating the wings almost horizontally with very high frequency, while the wing tip is

moving in a "figure 8" pattern. This is called "insect stroke" and is shown in Fig 1-2(B).

Take-off flight. This is the portion of the flight in which birds need to acquire

enough speed to utilize the aerodynamic forces for a normal flight. Small birds can use

their legs to jump in the initial stages of take-off, whereas large heavy birds take long

runs against the wind to reach take-off speed.

Landing. In this flight mode, birds need a large drag force to bleed their speed,

while maintaining the lift necessary to sustain their weight at low speed. The wing shape

during landing prevents flow separation at high angle of attack, and usually the tail is

used as a flap to increase lift.

1.2.3 Beating Frequency and Scaling

The mode and frequency of the beating motion depend on body size and shape,

flight mode, and flight regime. Therefore, different species use very different beating

motions. One thing they share is the ability to select the flapping mode and frequency that

lead to optimal power consumption of the respective flight mode for normal flight.

downstroke (power stroke)

upstroke




downstroke


upstroke (recovery stroke)
A) B)

Figure 1-2. Hovering flight. A) Avian stroke. B) Figure "8" pattern observed in hovering
flight of hummingbirds and most of the insects. Sketched from Azuma (1992).









The range of beating frequency and Reynolds number varies over a large range (as

seen in Table 1-1), where the body mass, beating frequency and Reynolds number are

presented for different species of birds and insects. In hovering flight, the maximum

velocity of the wing tip is given by the formula (Azuma, 1992):

U f -b- (1-1)

where f is the flapping frequency, b is the wing span, and /7 is the flapping

amplitude. The maximum Reynolds number at the tip is defined as (Azuma, 1992):

Re = U/uv = ((7fb/cluv) = (#fb2 /vAR)f, (1-2)

where is the mean chord defined by c = b/AR, and AR is the wing aspect ratio.

Then the reduced frequency can be defined as (Azuma, 1992):

k =rf c/U= 1/(AR 1) 1/AR (1-3)

The above equation states that a bird or an insect with a high aspect ratio wing flies

at a small reduced frequency. Generally, larger flying creatures fly by gliding or slow

beating, while smaller ones fly by strong beating at high frequency.

Two main factors limit the beating frequency in flapping flight. One is the

structural limitation of the bones (for birds), the supporting "veins" (for insects), or the

muscles, wherein a higher load can lead to structural failures. Another constraint is the

power available from the muscles. These limitations create the upper and lower bound for

the wing beat frequency. The power consumed to overcome the inertial torque of the

wing is proportional to the product of the moment of inertia and the cube of the beating

frequency as shown below (Azuma, 1992):

P I. f3 m5/3 f3 (1-4)









Table 1-1. Variety of beating frequency and Reynolds number for birds and insects.
Species Mass Beating Reynolds
m [kg] frequency, [Hz] number, Re
Chalcid wasp .
Chalcid wasp 2.5x10-8 370 2.0x101
(Encarsia formosa)
Fruit fly 2.0x10-6 240 2.0x102
(Drosophila virlis)
Fly 1.2x10-' 190 6.4x102
(Musca domestic)
Mosquito 3.5x10-6 320 1.3x103
(Aedes nearcticus)
Honeybee 7.8x10-5 250 2.3x103
(Apis mellifica L.)
Bumblebee
Bumbutlebee 8.8x10-4 156 4.0x103
(Bombus terrestris)
Dragonfly 7.9x10-4 29 4.9x103
(Anax parthenope)
Locust
Locust .2.0x10-3 20 1.0x104
(Schistocerca gregaria L.)
Hummingbird 2 4
Hummingbird 2.0x10-2 15 1.5x104
(Patagona gigas)
Sparrow 3.0x10-2 13 1.0x105
Pigeon 3.5x10-1 6 2.0x105
Stork 3.5 2 4.0x105
Data collected from Azuma (1992), Osbourne (1951) and Jensen (1956).

Azuma (1992) states that the maximum power generated by a muscle during one

cycle of flapping can be considered proportional to the product of the mass and the

frequency:

P a m.f (1-5)

Equating the maximum available power and the consumed power, the upper limit

of the wingbeat frequency is given by:

Pn oc P => m5/3 f3 o m f => fmax C m1/3 (1-6)

This same relation was specified by Hill (1950), who was one of the first to observe

the direct proportionality between the mechanical power produced by a particular flight









muscle and the contraction frequency. This has made flapping frequency an important

parameter when trying to describe the theory behind flapping wings.

To find the lower limit of wingbeat frequency, one can consider the force

equilibrium in cruising flight as the flight mode with minimum power consumption, and

consequently, with the lowest flapping frequency. In the vertical direction, the lift force

cancels the weight and the following proportionalities can be written:

L = W >pU2SC = mgg = U2S oc m => Uoc ~S => U oc 1/6 (1-7)


where the geometric similarity for the wing area S was used (S c m2/3) (Azuma,

1992).

In the horizontal direction, the thrust and drag are in equilibrium and the following

proportionalities are obtained:

1 1/3
T=D=> P/U= pU2SC => S(bf)3/UoU2S=> f ocUb => f ocUm1, (1-8)
2

where PA is the aerodynamic power and the geometric similarity for wingspan b was used

(b oc m13) (Azuma, 1992). Combining equations (1-7) and (1-8), the lower bound for

wingbeat frequency is given as:

f OcUm 13 o m1/6 (1-9)

The same expression for lower bound of frequency was obtained by Shyy et al.

(1999), starting with the assumption that the lower frequency is experienced by most of

the birds in slow forward flight or hovering. Pennycuick (1996), using dimensional

analysis, gave the following relation between the beating frequency f wing span b, and

the wing area S for birds in cruising flight:

f o m3 /b 23/24 1/3 3/8gl/2 (1-10)










where p is the air density and g is the gravitational acceleration. Introducing geometric

similarity expressions for the wing span and wing area (b oc m/3 S ocm2/3) (Azuma,

1992), equation (1-10) leads to: f oc m1/6, which is consistent with the result given in


equation (1-9).

Based on statistical data from birds and insects, Greenewalt (1962) showed that the

wingbeat frequency, f varies with the wing length lw (~ b/2) according to the relation:

f -1 = constant. For insects, the stroke frequency is widely scattered when plotted


against the body mass but is confined between the m-1/2 and m-1/6 line, as can be seen in

Fig. 1-3.


4Wing beat frequency vs body mass for different insect species
10
,- ,,,,ai different insect species
-I- ....... .. .- f = -1/2
--4-41+ 1 Ts 1 44f ff1-- 4
3I 1 1 1 1 11-ll I I I 1/3

1 I l 1 1 1 1 I I I Ill f f I
N 110 1















10 10 10 10 10 10-3 10
l Lij ^ J _~ 10- fii l6_II11__Illl_..









Figure 1-3. Wing beat frequency versus the body mass for different insect species. (Data










collected from Azuma, 1992).
I ----I 4IT -- I--1--I4I1 -- -- I7-44 -I -- F F1T IT -I-- T IT I -- I4 --I-- T I--I I-
Il 2 I I I Illi ll I I I Il l l l l I I I I l l l l
C81 -FF-Ii Ti T ii _L I I I I LI I III
S I- 1 1 T1 1 1 F T F i T T T I I F 4 III -1 F TI IT I II I II I I I IIIII I I I IIIII
l< 1 2 1 0 2 oi lI ,.II I I II I I I I I I II I I I I III lI l ll l I I I I I I II


I I I-II -4II I-I4 4III I4II I I I I I I -I TITI II III




108 10- 7 10T6 105 10- 4 10 3 102


Figure 1-3. Wing beat frequency versus the body mass for different insect species. (Data









flapping flight, P, and available power, P,, increases with the mass of the flying animal,









and therefore confirms the fact that the size of flying birds and insects is inevitably

limited.

A review of flapping wing kinematics and aerodynamic models for analyzing the

lift, drag and power at low Reynolds numbers, as well as the role of flexible structures

capable of adjusting to the free stream variations can be found in Shyy et al. (1999).

1.3 Insect Flight

The main objective of this work is to develop a better understanding of the flapping

wings aerodynamics at low Reynolds numbers, for possible applications to micro air

vehicle design. Therefore, the main focus of this study is on insect flapping wings. First,

the main differences between birds and insects flight regimes are presented. Second, the

main efforts in experimental research of the insect flight mechanisms are briefly

reviewed. The immense progress in computational methods and computer power allows

us to perform simulations of the flow over an insect flapping wing, and get more insight

into the unsteady flight mechanics of insect flight. These efforts are presented next,

followed by a description of the main mechanisms responsible for the aerodynamic force

generation.

1.3.1 General Characteristics

Based on their sizes, flying insects operate over a range of Reynolds numbers from

approximately 10 to 103. High beat frequencies and small sizes (0.25 mm to 0.1 m) result

in low wing loadings (less than 0.1 N/m2) and high reduced frequencies (more than 0.1).

Even if the inertial forces are at least one order of magnitude larger than viscous forces,

viscous effects have an important role in structuring the flow, and thus cannot be ignored.

Birds and insects use the aerodynamic forces differently. Flying at high Reynolds

numbers, birds only use the lift force, while reducing drag. Insects use both lift and drag,









which are strongly coupled by unsteady effects generated by the low Reynolds number,

high flapping frequency, and relatively low wing aspect ratio.

Generally, insects show remarkable flight agility and maneuverability, hover for

long periods, fly in any direction (even backwards), and are capable of landing upside

down on leaves or flowers (Alexander, 1986; Azuma, 1988).

For an airfoil, as the Reynolds number decreases, the maximum CL CD ratio

decreases and aerodynamic performances deteriorate. For birds, the feathers overlap and

form a profile that resembles a smooth "technical" airfoil. Insect airfoils have rough

surfaces such as well-defined cross-sectional corrugations (dragonfly) or scales on the

wing surface (butterfly and moth). Nachtigall (1976) reported that the removal of the

scales from a moth flying at a Reynolds number of 3x103 reduced lift by an average of

15% compared with the scaled wing. Kesel (2000) compared the aerodynamic

characteristics of dragonfly wing sections with technical airfoils and flat plates. In her

experimental study, she concluded that corrugated airfoils, such as those seen in

dragonflies (Fig. 1-4), have very low drag coefficients closely resembling those of flat

plates, while the lift coefficient was much higher than those of flat plates. The

measurements were done at two Reynolds numbers: 7.8x103 and 1x104. The cross-

sectional configurations vary greatly along the longitudinal axis of the wing, producing

different aerodynamic characteristics.

1.3.2 Experimental Efforts in Studying the Insect Flight Mechanisms

Because of their high wing beat frequencies and small size, it is quite difficult to

measure the wing motions of free flying insects. Marey (1868) was the first to use slow-

motion photography to track tethered insects' wingtip motion. More recent methods have









employed high-speed videography (Willmott and Ellington, 1997b), which offers greater

light sensitivity. Measuring the aerodynamic forces history on insects is a very

challenging task. At best, the forces have been measured on the body of a tethered insect,

making it very difficult to separate the inertial forces from the aerodynamic forces

(Cloupeau et al., 1979; Buckholz, 1981; Wilkin and Williams, 1993).



















Figure 1-4. Cross-sectional corrugations of dragonfly wings. Left wing is the hind wing
and to the right is the fore wing. Picture taken by Dragos Viieru.

Experiments on tethered insect flight give an important insight into the flow

structures in the wake, but due to their dimensions the visualization results of the flow

near the flapping wings are not clear enough, and in addition, tethering can alter the wing

motion as compared with free flight conditions. Also, it is very difficult to measure the

forces and moments generated by the insect wings in all degrees of freedom, and with

enough temporal resolution required to link the wake structures with the aerodynamic

force generation. To overcome these difficulties, many researchers have developed

mechanical models of insect wings, matching the Reynolds number and reduced

frequency parameter (body velocity/wing velocity) to that of an actual insect. These









robotic devices have provided new insights into the vortex structures and aerodynamic

forces of insect flight. A summary of the mechanical and robotic models is presented in

Table 1-2.

1.3.3 Computational Efforts in Studying the Insect Flight Mechanisms

The recent advances in computational methods and increasing computer power

offer a new approach in solving the mechanism of insect flight. The computational

methods applied to solving the fluid flow equations are more rigorous than the simplified

analytical solutions, but still require empirical data for validation of wing kinematics.

Nevertheless, based on more relevant input data from experiments, the computational

fluid dynamics (CFD) models offer a better insight in the unsteady mechanisms of

flapping flight.

One of the first computational studies modeling the flight of a tethered hawkmoth

was done by Smith (1996). In his investigation, he used an unsteady panel method for

solving the governing equations and accounted for the wing flexibility using a finite

element model. The aerodynamic forces at the center of the aerodynamic panels were

transferred to the structural control points using a load transformation matrix computed

using a finite surface spline method (Appa, 1989). The simulation results showed the

effects of both the aerodynamic and inertial forces on flapping, flexible wings. The study

also established the importance of including the wake in analysis of the unsteady

aerodynamics of flexible wings undergoing arbitrary large-scale motion.

Liu and Kawachi (1998) were the first to attempt a full unsteady Navier-Stokes

simulation of the flow over a hawkmoth wing on a structured grid. To avoid the artificial

mass fluxes due to grid movement, the "Geometric Conservation Law" (Thomas and










Lombard, 1979) was implemented. The kinematic model of the flapping wing was based

on the kinematic analysis of a hovering hawkmoth by Willmott and Ellington (1997a).

Table 1-2. Experimental robotic and mechanic devices for insect flapping wing study.
Author, Year Wing model Main features
based on:
Bennett (1970) Beetle 3-D mechanical model
(Melolontha flow visualization
vulgaris) only wing incidence can modified
do not allow force measurements
Spedding and Generic 2-D mechanical model for fling mechanism
Maxworthy (1986) only wing rotation allowed
allow measurements of the aerodynamic forces using a force
transducer
Saharon and Luttgers Dragonfly 3-D mechanical model
(1987), (1988) and flow visualization
Savage et al. (1979)
Ellington et al. (1996) Hawkmoth 3-D robotic model
good flow visualization (model wing =-10 times insect wing)
computer-controlled wing kinematic
do not allow force measurements
Dickinson et al. (1999) Fruit fly 3-D robotic model
and Lehmann (2000) good flow visualization (scaled-up wing)
computer-controlled wing kinematic
allow measurements of shear forces using a 2-D force
transducer

To facilitate the resolution of viscous flow at the leading edge, trailing edge and

wing tip, an 0-0 type grid topology was used. This study confirmed the smoke streak

patterns observed in both real and dynamically scaled model insect flight (Ellington et al.,

1996), and accurately predicted the complex vortex structure and the importance of the

span-wise flow in stabilizing the spiral leading-edge vortex as a lift-enhancement

mechanism, as one can see in Fig. 1-5.

Ramamurti and Sandberg (2002) studied the unsteady flow over a rigid fruit fly

wing based on a dynamically scaled model using a finite element solver and a re-meshing

technique. A detailed description of the mesh movement algorithms is given in

Ramamurti et al. (1994). The effect of phasing between the translational and rotational

motions was studied.


































Figure 1-5. Leading edge vortex on a hawkmoth wing. A) Numerical result. B) Smoke-
visualized flow during downstroke. Reprinted from Liu et al. (1998) with
permission from The Company of Biologists Ltd.

It was observed that the highest peak in the thrust force is obtained when the wing

rotation is advanced with respect to the stroke reversal, followed by the case in which the

wing rotation is symmetrical with respect to the stroke reversal. The lowest thrust is

obtained when the wing rotation is delayed. The computed thrust and drag roughly

matched the experimental results measured on a dynamically scaled model (Dickinson et

al., 1999), and showed the importance of the rotational mechanisms in accurately

describing the force time histories.

Another analysis on the unsteady aerodynamics of a model fruit fly wing, including

the energy consumption aspect, was conducted by Sun and Tang (2002a, 2002b) using an

implicit algorithm that is second-order accurate in time and space on a structured grid.









The wing movement was handled using a time dependent coordinate transformation

between an inertial reference frame and the non-inertial, body fixed frame of reference.

Three mechanisms were observed to account for the large lift: the rapid acceleration of

the wing at the beginning of a stroke, the delayed stall during the stroke, and the fast

pitching-up rotation of the wing near the end of the stroke (advanced rotation). The

computed unsteady lift coefficients were smaller than those measured on a robotic wing

by Dickinson et al. (1999). This was partly explained by using a wing with a smaller

aspect ratio than that of the robotic wing. However, the overall computed lift coefficient

pattern was very similar to the experimental results.

Recently, Gilmanov and Sotiropoulos (2005) performed simulations of flow over

the fruit fly (Drosophila). The wing kinematics and flow parameters were those

considered in the experimental work of Birch and Dickinson (2003). Their results showed

a reasonable agreement between the computed force history and the experimental values

(Birch and Dickinson, 2003). They numerically solved the three-dimensional, unsteady,

incompressible Navier-Stokes equations in a Cartesian domain using the immersed

boundary method, which is capable of handling complex geometries moving with

prescribed kinematics. The surface of the body moving in a fixed Cartesian domain is

discretized using an unstructured triangular mesh. The nodes of the surface mesh

constitute a set of Lagrangian control points used to track the movement of the body. At

each time step the influence of the body on the flow is accounted for by applying

boundary conditions at Cartesian grid nodes located in the exterior, but in the immediate

vicinity of the body by reconstructing the solution along the local normal to the body

surface. To discretize the governing equations, second-order finite difference formulas









are used on a hybrid staggered/non-staggered grid layout. The discrete equations are

integrated in time using a second-order, dual-time-stepping, artificial compressibility

iteration scheme.

The aerodynamic force generation and mechanical power requirements of a

dragonfly in hovering flight were studied by Sun and Lan (2004). Two O-H type grids

were generated around the forewing and the hindwing, and they moved in a fixed

Cartesian grid. The wing kinematics were based on experimental data of a hovering

dragonfly (Norberg, 1975; Reavis and Luttges, 1988). The computational results showed

two large vertical force peaks in one flapping cycle generated by a new vortex ring

containing downward momentum forming on the forewing or the hindwing during each

downstroke. The interaction between the forewing and hindwing was not found to be

very strong, and even detrimental to the vertical force generation. The study also showed

that the dragonfly uses drag as a major source of vertical force (approximately 65% of the

total vertical force is contributed by the drag) during hovering. Similar results were

obtained by Mittal et al. (2002) in their two-dimensional simulation of a hovering

dragonfly. The solver employed a second-order accurate central difference scheme for

spatial discretization and a mixed explicit-implicit fractional step scheme for time

marching (Udayakumar et al., 1999; Ye et al., 1999; Udayakumar et al., 2001). Different

values of phase lag between the forewing and hindwing were used. The results showed

that, for the set of simulations studied, the overall efficiency of the paired wing system

was lower than a single wing for all values of phase lag. For the single wing, the

simulation indicated that the highest thrust efficiency was accompanied by the formation

of an inverse Karman vortex street.










Despite the importance of 3-D effects, comparison of experiments and

computations in 2-D have also provided important insight. Wang (2000a, 2000b)

demonstrated that the force dynamics of 2-D wings might be sufficient to explain some of

the enhanced lift coefficients measured in insects.

A summary of the main numerical studies of the insect flapping flight and

governing equations solved are presented in Table 1-3.

Table 1-3. Principal numerical studies of insect flapping flight
Author, Year Fluid governing equations Method Moving boundary tracking
Smith, (1996) 3-D potential flow Panel method Time depended boundary
conditions
Liu and Kawachi, 3-D Navier-Stokes in Artificial Grid re-meshing.
(1998) body-fitted coordinates compressibility Analytical based on the initial
(Chorin, 1968) grid and wing kinematics
Wang, (2000a, 2000b) 2-D Navier-Stokes for Finite difference Time depending boundary
vorticity in elliptic conditions
coordinates.
Mittal et al., (2002) 2-D Navier-Stokes Finite volume Immersed Boundary Method
equations on fixed Cartesian grid
Ramamurti and 3-D Navier-Stokes in Finite Element Grid re-meshing
Sandberg, (2002) ALE formulation Spring analogy and smoothing
(Arbitrary Lagrangian on unstructured grid
Eulerian)
Sun and Tang, (2002a, 3-D Navier-Stokes in Artificial Time dependent coordinate
2002b) body-fitted coordinates, compressibility transformation
(Rogers et al.,
1991; Rogers and
Kwak, 1990)
Sun and Lan, (21 114) 3-D Navier-Stokes in Artificial Overset moving grid
body-fitted coordinates, compressibility
(Rogers et al.,
1991; Rogers and
Kwak, 1990)
Gilmanov and 3-D Navier-Stokes Finite difference Immersed Boundary Method
Sotiropoulos, (2005) on fixed Cartesian grid
Mittal et al., (2006) 3-D Navier-Stokes Finite difference Immersed Boundary Method
on fixed Cartesian grid


1.4 Unsteady mechanisms identified to enhance lift in insects

1.4.1 Analytical and quasi-steady models of insect flight

Early models of insect flight analyzed the far-field wakes and could not be used to

calculate instantaneous forces on airfoils, but could characterize the average forces. The









most notable analytical models are the "vortex models" of Ellington (1978, 1980, 1984e)

and Rayner (1979a, 1979b) that are applicable for hovering and forward flight. These

vortex theories provide a method to compute the mean lift and induced power for a given

circulation profile during the wingbeat. They consider the flapping wings as an actuator

disc (pulsed actuator disc) that pushes air downwards at a certain rate. The rate of change

of momentum flux within the downward jet must be equal to the insect weight, and

therefore the circulation in the wake required to maintain the force balance can be

calculated, leading to an estimation of the mean lift force. However, the momentum jet

theory touches the insect flight on a superficial level, ignoring variations in kinematics

and lift generation. A detailed description of these theories is presented in Rayner (1979a,

1979b) and Ellington (1984e).

Azuma (1985) and Azuma and Watanabe (1988) used the "local circulation

method" to study dragonfly and damselfly flight. This method incorporates the temporal

and spatial (span-wise) changes in induced velocity and corrects the circulation due to

wake. A more recent analytical model for flapping flight, assuming rigid wings, is

presented in Zbikowski (2002).

To predict aerodynamic forces from the motion of the wing, scientists developed

simplified models based on the quasi-steady approximation. The main assumption is that

the instantaneous aerodynamic forces on the flapping wing are equal to those generated

during steady motion of the wing at identical instantaneous velocities and angles of attack

(Ellington, 1984a). Therefore, the time dependence of the aerodynamic forces arises from

the kinematics' time dependence, and not from changes in the fluid flow. Also, the lift

and drag coefficients are time invariant and derived typically from experiments in which









model wings are tested over the range of angles of attack and flow velocities encountered

during flapping motion. The theory was used with some success by Osborne (1951) and

Jensen (1956) in their study of the desert locust's forward flight. However, the existing

quasi-steady theory predicts lower lift coefficients than the required average lift

coefficients for hovering as argued by Ellington (1984a, 1984b, 1984c, 1984d, 1984e) in

a comprehensive review of the insect flight literature and based on a wide range of data

available at that time. Direct force measurements in flying insects also validated his

findings, indicating that insects do not use the steady-state principle (Ellington, 1984a,

1984b, 1984c, 1984d, 1984e). Based on experimental, computational and analytical

models, some important unsteady mechanisms that play an important role in generating

aerodynamic forces in insect flight can be identified, and they are described below.

1.4.2 Wagner Effect

The Wagner effect is manifested when an inclined wing starts impulsively from

rest. In this case the circulation around the wing does not immediately attain its steady

state value. This mechanism that delays the development of the maximum circulation was

proposed by Wagner (1925) and studied experimentally by Walker (1931). Based on

Walker's results, the steady circulation is established after approximately 4 to 5 chord

lengths. Ellington (1984d) defined the ratio of the mean distance traveled by the wing

during a half-stroke to the mean chord as a measure of the number of chord lengths

traveled. For the most insects this value is between 1.9 and 3.8; this fact indicates that

steady circulation is never reached (Ellington, 1984a). It seems that, unlike other

unsteady mechanisms, the Wagner effect would attenuate the unsteady forces below

levels predicted by quasi-steady models. Studies done on 2-D wings by Dickinson and

Gotz (1993) indicate that the Wagner effect is not very strong in the range of the









Reynolds numbers experienced in insect flight, and most recent models of insect flight

neglected the Wagner effect.

1.4.3 Leading-edge Vortex and Delayed Stall

For thin airfoils traveling at high angle of attack the local viscous forces within the

fluid near the leading edge are smaller than the pressure forces generated by the high

fluid velocity, and the flow separates from the upper surface of the wing. The low

pressure region behind the leading edge forces the flow to curl back and reattach on the

upper wing surface behind the leading edge, leading to the formation of a leading edge

vortex. In this case a suction force develops perpendicular to the wing chord,

consequently enhancing the lift as well as the drag, as shown in Fig. 1-6.

For a 2-D motion of the wing at high angle of attack, the leading edge vortex grows

in size until reattachment is no longer possible. The leading edge vortex sheds into the

wake, and as a result the lift drops and the wing stalls. But for several chord distances

before stall, the attached leading edge vortex generates very high lift coefficients

(delayed stall). At low Reynolds numbers, as the leading edge vortex is shed, another

vortex is formed at the trailing edge that grows until it is also shed into the wake and a

new leading edge vortex is forming.

This dynamic process creates a wake of regularly spaced counter-rotating vortices

known as the Karman vortex street (Dickinson and Gotz, 1993). One of the first

visualizations of the leading edge vortex and the delayed stall mechanism in insect flight

was done by Maxworthy (1979) on the model of a flinging wing. Since most insects flap

their wings at high angle of attack, the leading edge vortex has an important role in lift

generation.









FResultant Fsucton + FNormal
Lift



Fsuction




Drag







Figure 1-6. Flow around a thin airfoil. The presence of the leading edge vortex generates
a suction force normal to the surface of the airfoil due to the low pressure
zone associated with the vortex core. This force enhances the resultant force
normal to the airfoil. The filled circle indicates the leading edge. Sketched
from Sane (2003).

In contrast with 2-D models, the leading edge vortex in 3-D flows was not shed

even after many chord lengths and no Karman vortex street is present, as observed by

Ellington et al. (1996) in his smoke visualization of the flow around both real and robotic

hawkmoth wing's at a Reynolds number in the range of 103. The mechanisms of 3-D

leading edge vortex stability are still unclear. Ellington and co-workers observed a steady

span-wise flow from the wing's hinge to approximately three-quarters of the span

towards the tip where the conical leading edge vortex detaches from the wing surface.

This phenomenon was also observed by Liu and Kawachi (1998) in their computational

study of the hawkmoth wing. The axial flow role in stabilizing the leading edge vortex

was also observed in the case of back-sweep delta wings. However, the importance of

axial flow in the stability of the leading edge vortex was questioned by experiments using

a 3-D model wing based on the fruit fly flapping at a Reynolds number in the order of

100 (Birch and Dickinson, 2001). The measurements showed that the span-wise flow









within the leading edge vortex core was quite small (2-5% of the average tip velocity),

while a large axial flow on the wing's upper surface behind the leading edge vortex was

observed, that rolls into a large tip vortex.

Recently, in an extended study of the effect of Reynolds number on the leading-

edge vortex at mid-downstroke in hovering, Liu et al. (2005) found that the leading-edge

vortices show significant variations in their structures as the Reynolds number changes.

Figure 1-7 shows the streamline structures at three Reynolds numbers, corresponding to a

hawkmoth in Fig. 1-7(A), a fruit fly in Fig. 1-7(B), and a thrips in Fig. 1-7(C). At a

Reynolds number of 6,000 (corresponding to a moth) as observed in the previous studies

(Ellington, 1996; Liu and Kawachi, 1998), an intense, conical leading-edge vortex is

observed on the paired wings with a steady span-wise flow at the vortex core, breaking

down at approximately three-quarters of the span towards the tip. At a Reynolds number

of 120, corresponding to a fly (Fig. 1-7(B)), the break-down apparently disappears, and a

leading-edge vortex is found to be connected to the tip vortex, forming a longer, leading-

edge tip vortex. The vortex shows strong stability, while the span-wise flow at the vortex

core turns to become much weaker compared with those at high Reynolds numbers,

which is in agreement with the measurements of Birch and Dickinson (2001). Further

reducing the Reynolds number down to 10 for the thrips case, a vortex ring connecting

the leading-edge vortex, the tip vortex and the trailing vortex is observed (Fig. 1-7(C)).

Their 3D structure looks more like a cylinder than the conical form seen at high Reynolds

numbers. Inspecting the momentum equation, one can see that the pressure-gradient, the

centrifugal force, and the coriolis forces are all together responsible for the leading-edge









vortex stability. Their roles at different Reynolds numbers should be studied to help shed

light on these findings.











Figure 1-7. Leading edge vortex (LEV) structures at different Reynolds numbers. A)
Reynolds number = 6000. B) Reynolds number = 120. C) Reynolds number =
10. Adopted from Liu et al. (2005).

In another effort, Sane and Dickinson (2002), by placing a fence on the upper

surface of a robotic model wing, showed that the leading edge vortex (LEV) still exists

on the wing. Based on this observation, they suggested that the spanwise flow should not

be responsible for the LEV on the upper wing surface.

1.4.4 Rotational Forces

When a wing rotates about a span-wise axis and translates at the same time, it

generates a rotational circulation in the fluid to counteract the effects of the rotation. The

additional circulation is proportional to the angular velocity of the rotation and generates

rotational forces, which, depending on the direction of the rotation, can be added or

subtracted from the net force due to translation. Insects rotate their wings around the

span-wise axis near the end of every stroke in order to obtain a positive angle of attack

for the next half-stroke and thus, the forces generated during rotation are very important.

Direct force measurements on tethered fruit flies have shown that maximum forces are

produced during wing rotation, suggesting that the rotational forces generated during

stroke reversal represent a significant part of the total lift (Dickinson and Gotz, 1996).









Similar rotational force peaks were observed in CFD simulations by Sun and Tang

(2002a). Both CFD simulations and experimental force measurements indicate that the

forces during flight are sensitive to the phase of wing rotation. Advancement in wing

rotation at the end of a half-stroke generates the maximum mean lift coefficient, while

delaying the rotation reduces the mean lift coefficient by almost 42% (Dickinson et al.,

1999). Dickinson et al. (1993) suggested that rotation timing is used by insects during

steering maneuvers; the insect advances the rotation on the wing outside the turn, and

delays it on the wing inside the turn.

1.4.5 Wing-wake Interaction

In contrast with wings of aircraft that move through still air, the wings of hovering

insects intercept the wake created by the wing's own movement. The interaction of insect

wings with the shed vortex of prior strokes was observed with both force measurements

and flow visualization on a robotic model based on a fruit fly (Dickinson et al., 1999).

During stroke reversal both leading and trailing edge vortices are shed, and consequently

a strong velocity field is generated between these two vortices. As the wing reverses

direction, it encounters the enhanced velocity field that transfers momentum to the wing,

thereby enhancing the maximum force during stroke reversal. This mechanism is

illustrated in Fig. 1-8. The magnitude and relative strengths of shed vortices are strongly

dependent on the wing kinematics immediately before and after stroke reversal.

The model wing produces positive wake capturing lift when the wing rotation

precedes stroke reversal. Similar force peaks were observed by Sun and Tang (2002a) in

their CFD simulations of a fruit fly based on kinematics similar to those in Dickinson et

al. (1999), but they attribute these peaks to the fast pitching-up rotation of the wing, and

not to the wing-wake interaction. They also provided some experimental cases to prove









their conclusion. Later, Walker (2002), using a different CFD method, reached the same

conclusion.

A
D





B E Transferred
\ momentum

O Downwas
Local flow

C O F






Figure 1-8. Momentum transfer due to wing-wake interaction. A) Wing is in steady
translation phase. B) Leading and trailing edge vortices are generated as the
wing rotates around a spanwise axis. C) Velocity field induced by the shed
vortices. D) Wing rotates and starts the reverse stroke. E) Wing reverses the
stroke and encounters the induced velocity field. Fluid momentum is
transferred to the wing that generates a peak in the aerodynamic. F) Wing is in
reversed stroke steady translation. The filled circle indicates the leading edge.
Sketched from Sane (2003).

The wake capturing mechanism may explain the insect's ability to extract energy

from its own wake, recovering energy from air lost during the previous stroke and

therefore improving the overall force production efficiency. Srygley et al. (2002)

observed, during a smoke visualization experiment, the extraction of kinetic energy from

the wake behind freely flying butterflies during take-off. However, the contribution of the

wing-wake interaction mechanism to the peak force generation during wing rotation is

still not clear, and requires further investigation.









1.4.6 Clap and Fling Mechanism

The clap-and-fling mechanism is the physical interaction of the left and right wing

during dorsal stroke reversal. This mechanism was proposed by Weis-Fough (1973) to

explain the lift generation in the chalcid wasp. Its schematic is presented in Fig. 1-9.



clap

a b c d

fling v


a b c d

Figure 1-9. Schematic of clap and fling mechanism. Based on Weis-Fough (1973).

This mechanism was also found in fruit flies (Gotz, 1987), butterflies (Brodsky,

1991), and locusts (Cooter and Baker, 1977). During the fling phase, preceding the down-

stroke, the fluid is pulled into the gap formed by the moving wings, generating strong

leading edge vortices while the development of the trailing edge vortices is stopped by

the contact of the wing's trailing edges. When the wings translate away from each other,

the low pressure region between them accelerates the surrounding fluid, providing a

build-up of circulation. A detailed theoretical analysis of the clap-and-fling mechanism

can be found in Lighthill (1973) and Sunada et al. (1993), and experimental observations

can be found in Bennett (1977), Maxworthy (1979), and Spedding & Maxworthy (1986).

The overall benefit of the clap and fling remains uncertain in many studies of insect wing

kinematics.

Recently, the unsteady mechanisms of lift enhancement in insect flight were

reviewed in Lehmann (2004). A review of the current models as well as of the unsteady









mechanisms of insect flight is presented in Sane (2003). Recent findings in the unsteady

aerodynamics of insect flight at intermediate Reynolds numbers (10-104), with focus on

two-dimensional unsteady mechanisms in uniform and accelerated motions, forward and

hovering flight can be found in Wang (2005). A comprehensive survey of the research of

flapping wing propulsors and their application to man-made vehicles has been presented

by Rozhdestvensky and Ryzhov (2003).

1.5 Objective and Outline

The main objective of this work is to investigate the aerodynamic force generation

mechanisms for fixed and flapping wings micro air vehicles at low Reynolds numbers. In

the light of the insect flight analysis presented above, it is clear that a better

understanding of the unsteady mechanisms for lift generation is very important if we

want to explore the development of micro and nano air vehicles. Therefore, in this work,

the focus is to better understand biological lift generation by flapping wings at low

Reynolds number, as well as to point out some of the characteristics of small micro air

vehicles with fixed rigid or flexible wings.

In this chapter, the flapping wing in nature was discussed with a focus on insect

flight. Main experimental and computational research efforts were presented, followed by

a presentation of the main unsteady mechanisms identified so far that are responsible for

lift generation in insects.

In chapter 2, the numerical framework for solving fluid-structure interaction

problems with moving boundaries is presented. The fluid's governing equations are

presented as well as the numerical implementation of those equations in a pressure-based

fluid solver. The finite element structural solver and the governing equations for solids

are also presented. The fluid-structure interaction problem leads to deformation in the









computational grid and therefore algorithms to automatically take into account the grid

movement are described. Since the fluid and structure solvers require different grid types,

an interpolation technique between the structured CFD grid and the unstructured FEM

(finite element method) mesh is mentioned. Finally, a generic algorithm for solving the

fluid-structure interaction problem is presented.

In chapter 3, the flow over fixed rigid and flexible wing MAVs is analyzed. The

low aspect-ratio wings used in the University of Florida MAV design generates strong tip

vortices that have an important impact on the wings' aerodynamic performances. The

wing tip vortex effect is studied along with means to reduce the induced drag and

increase the rigid wing efficiency. Lastly, the aerodynamic performances of a flexible

wing are compared with those of a similar rigid wing.

The two-dimensional hovering of an elliptic airfoil is investigated in chapter 4. The

results are compared with other computational work to assess the capabilities of the fluid

solver to handle this type of problem. Also, basic unsteady mechanisms for lift generation

are explored, such as delayed-stall and wake-capturing mechanisms.

Finally we conclude with the summary and the directions for the future work in

chapter 5.














CHAPTER 2
COMPUTATIONAL METHODS FOR COUPLED FLUID-STRUCTURE
INTERACTIONS

2.1 Overview

This study investigates the aerodynamics of fixed-wing MAVs with both rigid and

flexible wings, as well as the insect flapping flight under low Reynolds number

conditions. During flight, the flexible wing changes its shape as a response to the external

forces and, at the same time, its shape variation affects the pattern of its surrounding fluid

flow. Under steady freestream the membrane wing exhibits high frequency vibration,

observed both experimentally (Waszak et al., 2001) and numerically (Lian et al., 2003c;

Lian and Shyy, 2005). Also, the insect flapping flight is dominated by unsteady

mechanisms such as vortex shedding, dynamic stall, and wing-wake interaction. It is

desirable to have a computational capability to accurately capture the transient complex

behavior of fluid flow around flexible and flapping wings.

Numerous approaches have been proposed for solving unsteady flows. Pulliam

(1993) incorporated subiterations to improve the time accuracy of conventional implicit

schemes. Jameson (1991) and Rumsey et al. (1996) developed a dual time stepping

method used in the multigrid framework to improve the efficiency and accuracy of

traditional factored schemes. Issa (1985), and Oliveira and Issa (2001) proposed the

Pressure-Implicit Splitting of Operators (PISO) method to improve convergence. For a

review of these methods and some other frequently adopted approaches, see Shyy and

Mittal (1998).









In this work, the three-dimensional Navier-Stokes equations for incompressible

flow written in curvilinear coordinates are solved. Two popular methods to solve the

Navier-Stokes equations for incompressible flow are employed. One is the SIMPLEC

method (van Doormaal and Raithby, 1985) that is an enhancement of the SIMPLE

algorithm (Patankar and Spalding, 1972; Shyy, 1994), and the other one is the PISO

method (Issa, 1985). Both methods belong to the pressure-based family (Shyy, 1994) that

uses an artificially derived pressure or pressure correction equation. In this chapter, the

fluid solvers adopted for moving boundary problems are briefly introduced, as well as the

importance of the boundary conditions and the geometric conservation law (Thomas and

Lombard, 1979). Different moving grid techniques used to automatically re-mesh the grid

are investigated. Finally, the membrane structural solver and the coupling between the

fluid solver and structural solver are discussed.

2.2 Governing Equations for Fluid Flows

The three-dimensional Navier-Stokes equations for incompressible flow in

Cartesian coordinates can be written, using indicial notation, as follows:

ai+ a ( )=O (2-1)
at cx

Ca(u)+ I + u,) (2-2)
at ax ax, ax

where x, is the position vector, t is time, p is density, u, is the velocity vector, p is

pressure, and r, is the viscous stress tensor. The constitutive relation between stress and

strain rate for a Newtonian fluid is used to link the components of the stress tensor to

velocity gradients:









lu, = +Uj 2 u1 (2-3)
Oxi cx 3 Ox,

where / is the molecular viscosity.

For arbitrary shaped geometries, the Navier-Stokes equations are transformed into

generalized curvilinear coordinates (, r, C), where = (x,y,z), = r(x,y,z) and

S= $(x,y,z). The transformation of the physical domain (x, y, z) to the computational

domain (q, r, C) is achieved by the following relations (Tannehill et al., 1997):

[x Ay z 1 12 f 13
x y z = f21 f22 f23 (2-4)
S y z _f31 f32 f33

where f, are the metrics terms and J is the determinant of the Jacobian

transformation matrix given by:


S(x, y, z)x
J = det | -(Y-) = det yL Yq Yc (2-5)
_z z z _

To solve for the Navier-Stokes equations in curvilinear coordinates, the finite-

volume formulation is adopted. In this approach, both Cartesian velocity (as primary

variables) and contravariant velocity components are employed. The contravariant

velocities are used to evaluate the flux at the cell faces and to enforce the mass continuity

in the pressure-correction equation. The expressions for the metrics, the determinant of

the Jacobian transformation matrix, and the fluxes at the cell faces for the three-

dimensional Navier-Stokes equations in the generalized body-fitted curvilinear

coordinate system ( r), are given in Thakur et al. (2002) and Shyy (1994). When the









governing equations are considered under a moving grid framework, the grid velocities

should be included in the flux computations as described in Shyy et al. (1996).

2.3 Fluid Flow Solvers

The SIMPLEC method (van Doormaal and Raithby, 1985) used in this study to

solve the governing equations is an enhanced variant of the SIMPLE algorithm initially

proposed in Cartesian coordinates by Patankar and Spalding (1972). The extension to

body-fitted curvilinear coordinates capable of handling complex geometry computation

was discussed by Shyy (1994). An attractive feature of these methods is that the solution

procedure is sequential in nature and can therefore easily accommodate a varying number

of equations depending on the physics of the problem involved. The implementation of

the method employs a control volume approach and uses a non-staggered arrangement for

the velocity components and the scalar variables (i.e. pressure). The Cartesian velocity

components are computed from the respective momentum equations. The cell fluxes and

pressure fields are corrected using a pressure correction equation, which is derived by

manipulating the continuity and momentum equations. The iterative correction procedure

leads to a divergence-free velocity field within a desired convergence tolerance, therefore

enforcing the pressure-velocity coupling. The SIMPLEC method uses coordinated under-

relaxations for the momentum and pressure corrections to improve the convergence speed

(intrinsically slow in the original SIMPLE method).

The other method used in this work is the PISO method. Contrary to the iterative

SIMPLE family of methods, the PISO method uses a non-iterative process. Here the

pressure-velocity coupling is handled by splitting the solution process into a series of

predictor and corrector steps. The PISO method is generally more efficient for transient

flow computations (Issa, 1985; Thakur et al., 2002).









To solve the Navier-Stokes equations, proper boundary conditions are required. At

the interface between the fluid and solid structure, the no-slip condition is applied, which

requires that the fluid velocity at the wing surface matches the surface velocity of the

structure. These velocity boundary conditions are provided by computing the

displacements of the solid structure.

2.4 Geometric Conservation Law

To solve the governing equations in body-fitted curvilinear coordinates, a

transformation matrix is used to facilitate the mapping of a physical flow region (x, y, z)

onto a computational domain (4, r, Q. The Jacobian transformation matrix is defined as

in equation (2-5). The determinant of the Jacobian transformation matrix represents the

volume of the element in the transformed coordinate. In moving grid problems, the

computational grid is changing with time and consequently the determinant of the

Jacobian matrix, J, needs to be updated. Special procedures are required to compute the

effective value of J at each time step; otherwise errors arise due to an inconsistent

numerical implementation that would lead to the generation of artificial mass sources. As

suggested by Thomas and Lombard (1979), in the process of updating the Jacobian (J)

the following conservation law (GCL) needs to be satisfied:

d fJddqd;= f (V. W )d dd (2-6)
dtV V

where Vis the volume bounded by the closed surface S, and W, is the local velocity of the

moving boundary surface S. Thomas and Lombard (1979) proposed an expression to

evaluate the Jacobian (J) from the continuity equation for a constant density, uniform

velocity field under a time dependent coordinate transformation, while maintaining the

geometric conservation law:










= j J + 0 J (2-7)
Ot 9^ t0 B7 Wt 8 OR

Integrating equation (2-7) with a first-order, implicit time integration scheme over

the same control volume used for mass conservation leads to a finite volume discrete

form of the above equation. This is used to update the Jacobian in a manner that respects

the basic requirement of the geometric conservation in the discrete form of the

conservation law when the grid is time dependent. Implications of the geometric

conservation law on moving boundary problems were discuss by Shyy et al. (1996).

2.5 Moving Grid Technique

2.5.1 Master-slave Concept

For moving boundary problems where a solid boundary (i.e. wing) moves inside a

computational domain based on known kinematics (i.e. rigid flapping wing) or as a

response of the structure to the flow around it (i.e. membrane wing), the grid needs to be

adjusted dynamically during computation. To facilitate this, a moving grid technique

needs to be employed. The actual process of generating a grid is a complicated task, so an

automatic and fast algorithm to frequently upgrade the grid is essential. It is desirable to

have an automatic re-meshing algorithm to ensure that the dynamically moving grid

retains the quality of the initial grid, and avoids problems such as crossover of the grid

lines, crossed cell faces, or negative volumes at block interfaces (in case of the multi-

block grids). The moving grid technique can also be used for shape optimization studies

(Lian et al., 2003b).

Several approaches have been suggested to treat grid redistributions for moving

grid computations. Schuster et al. (1990) used an algebraic shearing method in their

study. The displacement of the moving surface is redistributed along the grid line which









connects the moving surface to the outer boundary. This simple method gives good

results for modest displacements and single block grids. For multiblock grid

arrangements extensive user intervention is required. A robust method that can handle

large deformations is the spring analogy method, which was first introduced by Batina

(1989) for unstructured grids and later extended to structured grids by Robinson et al.

(1991). In this method all edges of a cell, as well as the diagonals, are replaced by linear

springs, each with the stiffness inversely proportional to a power p of the length of the

connecting edges. Using the power p, one can control the stiffness of the spring and

consequently control the amount of movement to avoid excessive mesh distortions. The

iterative process necessary to find the displacement of all the internal points increases the

computational time of this method, especially for large grids. The direct transfinite

interpolation method was introduced by Eriksson (1982), and can generate grids for

complex geometry. The method defines an interpolation function given known values on

constant planes and function derivatives in the out of surface direction on the boundaries.

The method is fast and efficient for structured grids, but the quality of the initial grid may

not be always preserved, especially far from the boundaries. Hartwich and Agrawal

(1997) propose the Master/Slave concept to expedite the grid regeneration process and

preserve the grid continuity at the multiblock grid interfaces.

In this study a moving grid technique proposed by Lian (2003a) and Lian et al.

(2003d) is used to re-mesh the multiblock structured grid for fluid-structure interaction

problems. For multiblock structured grids, for simplicity, CFD solvers often require

point-matched grid block interfaces. This method is based on the master-slave concept









and maintains a point-matched grid block interface while maintaining grid quality and

preventing potential grid cross-over.

When an object changes its shape, the master points, which are located at the

moving body surface, move first, and then affect the distribution of the off-body points.

One difficulty for a multiblock grid resides in the way in which the vertices of each block

are moved, if a point-to-point match between two abutting blocks is required without

overlap. For identical interfaces between two abutting blocks, such as the interface

between block 2 and block 3 in Fig. 2-1(A), the edge and interior points can be obtained

by a 3-stage interpolation once the corer vertices are determined. However, when the

abutting blocks do not have an identical interface, such as of block 1 and block 2 in Fig.

2-1(A), the interpolation can cause a discontinuity at the interface. To avoid creating

undesirable grid discontinuities, the off-body vertices of a grid block are linked to a

surface point and thus move in a similar way. Therefore, for each off-body vertex (slave

point), the nearest body surface point is defined as its master point. The distance between

the slave point and its master is given by:


r = J(xs m)2 +(s Y)2 +s Zm (2-8)

where the subscript s represents a slave point, and m a master point. Once a slave point

has its master point identified, the slave point moves according to the influence from its

master. A simple formula to define the displacement is shown here:

Xs = Xs +0 (Xm xm) (2-9)

The tilde (-) indicates the new position, and 0 is a decay function used to control

and scale the movement to avoid a lower surface crossing over its opposite face when the

displacement is large.









i 1, )


RUM-I"'C" "
=- i:i: I%
-' 7: -
.7 :i.i..iiii:.i


Figure 2-1. Illustration of the moving grid technique. A) Initial geometry. B) Side view of
the initial geometry. C) Modified geometry. D) Grid points distribution of the
new shape. E) Same test but with higher stiffness. Adopted from Lian et al.
(2003c).

The decay function is a Gaussian distribution function given by:


0= exp{- p min[500, dv/( + dm)]l}


(2-10)


where e is a small number necessary to avoid division by zero when dm = 0, and


dv= (x, -xm )2 +(y 2 +(z Z- )2,

dm =(xm -m)2 +(Y _m)2 +(Z -Zm )2.


(2-11)











The coefficient f controls the stiffness of the movement. A larger value of f causes

the block to behave more like a rigid body as one can see in Fig. 2-1(E). More details

about master and slave concept and its applications can be found in Hartwich and

Agrawal (1997), Lian (2003a), and Lian et al. (2003c).



-0 18

-02 -
-0 22

0.21 -



-0,25
-026 -
-027
-028 i
-045 -0425 -0.4 -0.375 -035
x


Figure 2-2. Master-slave algorithm performance for large rotations

The master-slave algorithm is highly automated, capable of handling large

translational deformations while maintaining overall grid quality and allowing large

translational displacements. The propagation distance of the moving wall perturbation is

controlled by the spring stiffness coefficient. However, in the current formulation, large

rotational deformations cannot be handled properly since no information is provided

about cell skewness as the perturbation is propagated along grid lines. As one can see in

Fig. 2-2, for large deformations, the grid quality deteriorates, leading to negative cell

volumes.


2.5.2 Grid Re-arrangement Algorithm with Weighting Functions

The flapping wing problem requires both large rotations and displacements of the

airfoil, and an algorithm that can handle this type of problem needs to be developed. The

algorithm proposed dissipates the "shear stress" created by the movement of the body






40


toward the outer boundaries. A weight function is computed for each grid node based on

a weighted average of the neighboring node positions. The main steps of the grid-

rearrangement algorithm are as follows:

i) Assign the boundary conditions. Assign weight function (yp) = 1 (rigid

rotation/translation) on the body surface and 0 on the outer boundary (no displacement)

ii) Start the iterative process (k=l)

iii) Compute the weight function value:

Nnbr
Y (0.
O, = minbr where Nnbr is the number of neighboring nodes (8 for 2D -

m=l

see Fig. 2-3), / is the distance between the (i,j) node and the neighboring nodes, andp is a

stiffness parameter (default = 2, but can vary).

iv) check convergence criteria: (o (o, < tolerance. If YES => new weight

function, if NO => k=k+1

v) increment the new node position: x" ,' = x, + dx (, (translation) =>


x",j = x,""' + dy dac. (- (rotation), where dx, dy = translation increment, da = rotation

increment.


7 56 5

8 ij 4


m=l 2 3

Figure 2-3. Current node neighbors position for grid re-arrangement algorithm.







41


vi) at block interfaces the displacement is interpolated as the average between

neighboring blocks values, in order to maintain the point-to-point interface.

The main advantage of this algorithm is its ability to handle large rotational and

translational deformations, while maintaining a good grid quality near the body; as one

can see in Fig. 2-4(A).



14
i 12
"10 "







12
-1 14
0 1 -15 -10 -5 0 10 I 1
A) x B) x


Figure 2-4. Grid quality for the grid re-arrangement algorithm. A) Close to the body. B)
Overall domain.

The disadvantages of this algorithm are the large skewness near the outer boundary

(Fig. 2-4(B)) imposed by the non-movement condition and the approximation of dx, dy,

da quantities, leading to poor grid quality after a long time. This problem is resolved for

periodic movements by replacing the deformed grid with the original one after each

cycle.

2.5.3 Grid Re-arrangement Algorithm Based on Spring Analogy Concept

To overcome the main disadvantage of the grid re-arrangement algorithm described

in section 2.5.2 and of the master-slave algorithm (section 2.5.1), a new algorithm based

on the spring-analogy concept was developed by Tang et al. (2006). This algorithm

propagates the perturbation due to body movement along a grid line as described in










section 2.5.1, but also takes into account the angular deformation weighted with the

original grid orthogonality information. This way, the cross-over encountered in the

original master-slave algorithm (Fig. 2-2) is eliminated.

Table 2-1. Moving boundaries techniques. Advantages and disadvantages.
No. Grid Advantages Disadvantages Observations
movement
algorithm
1. Master-slave Can handle large translational Cannot handle large Used for MAV
(sec. 2.5.1) displacements rotation deformations studies small
Displacement can be controlled vertical disp.
by varying the stiffness
2. Weight Can handle large translational Error accumulation after Used for 2-D
function and rotational displacements. long time flapping wing
(sec 2.5.2) Displacement can be controlled studies (large
by varying the stiffness rotations).
parameter (p)
3. Spring Can handle large translational No current generalization Used for 2-D
analogy and rotational displacements for arbitrary multiblock flapping wing
(sec. 2.5.3) grids. studies (large
rotations).

Extensive computations also proved that the grid quality is maintained in time. The

main advantages and disadvantages of the grid movement routines are summarized in

Table 2-1.

2.6 Membrane Structural Solver

The main characteristic of membrane structures is their large deformations

implying a nonlinear stress-strain relation. Green and Adkins (1960) were the first to

derive a general nonlinear theory for membranes. The membrane theory is based on two

main assumptions (Jenkins and Leonard, 1991): the membrane tension stiffness is much

larger than its bending stiffness and the ratio between the membrane thickness and the

radius of the curvature is small (therefore the strain-displacement is decoupled from the

curvature tensor). Two approaches exist to model the response of a membrane structure.

The first approach idealizes a plate or shell structure as a membrane away from the

boundaries and neglects the stress couples in the interior region. This model can be









described by Foppl-von Karaman theory (Jenkins and Leonard, 1991; Jenkins, 1996). The

second approach assumes that the structure cannot sustain stress couples anywhere. This

approach is used in the current work and was also adopted by Lian et al. (2002) in their

simulation of a membrane wing.

The studies of a two-dimensional membrane with fixed leading and trailing edges

were originally conducted by Thwaites (1961) and Nelsen (1963). In these studies, an

inextensible membrane surrounded by steady, two-dimensional, irrotational flow was

considered. Under these assumptions, and for small camber and incidence angle, the

problem can be linearized and expressed in an integral form known as the "sail equation"

(Newman, 1987). A comprehensive review of earlier works related to membrane wing

aerodynamics can be found in Newman (1987).

Unlike the two-dimensional case, the three-dimensional membrane model

introduces several complicating factors. In the two-dimensional "sail equation", the

tension is a single constant. In three-dimensional membranes, the tension is defined as a

biaxial tension along the lines of principal stress (Jackson and Christie, 1987). Also, the

geometric and material properties may vary along the spanwise direction and need to be

described in detail. Another complicating factor is that the membrane can be compressed

and wrinkled when one of the principal tensions vanishes.

A finite element analysis of the static equilibrium of an inflated membrane

undergoing large deformations was presented by Oden and Sato (1967). A review of the

earlier work on the dynamic analysis of membranes can be found in Jenkins and Leonard

(1991). An update of the state of the art models in membrane dynamics is presented by

Jenkins (1996). Verron et al. (2001) have studied, both numerically and experimentally,









the dynamic inflation of a rubber-like membrane. Ding et al. (2003) numerically studied

partially wrinkled membranes.

The three-dimensional membrane model used in this work to study the membrane

dynamics has been proposed by Lian et al. (2002). The model gives good results for

membrane dynamics with large deformations but cannot handle the wrinkle phenomenon

that occurs when the membrane is compressed (Lian, 2003a). The membrane material

considered obeys the hyperelastic Mooney-Rivilin model (Mooney, 1940). A brief review

of their membrane model is given next.

2.6.1 Mooney-Rivlin Model

The Mooney-Rivlin model is one of the most frequently employed hyperelastic

models because of its mathematical simplicity and relatively good accuracy for

reasonably large strains (less than 150%) (Mooney, 1940).

For an initially isotropic membrane, a strain energy function, W, can be defined as

(Green and Adkins, 1960):

W= W(I,2,I3) (2-12)

where Ii, 12, and 13 are the first, second, and third invariants of the Green deformation

tensor, C. For an incompressible material, when 13 = 1, the strain energy is a function of

I, and 12 only, and a linear expression can be written for the membrane strain energy:

W =c,(, 3)+ c2 2- 3) (2-13)

where cl and c2 are two material constants. A material that obeys equation (2-13) is

known as a Mooney-Rivlin material. For an initially isotropic, incompressible membrane,

the general stress-strain relation is written as:









S= -p-C +2[(c, +c21)I-cC] (2-14)

where S is the second Piola-Kirchhoff stress tensor, p is the hydrostatic pressure, and I is

the 3x3 identity matrix. The membrane stress field is essentially assumed to be two-

dimensional and therefore, the stress normal to the deformed membrane surface is

negligible with respect to the stress in the tangent plane (Oden and Sato, 1967). Under

this assumption, the deformation matrix C(t) and the stress matrix S(t) can be written as:

C,,(t) C12(t) 0 S 1(t) S12(t) 0
C(t= C12 ( C22 (t) 0 ; S(t)= S12(t) S22(t) 0
0 0 C33(t) 0 0 0. (2-15)

The hydrostatic pressure is determined by the condition that S33 = 0, and the formula is:

p=2(c +C2I1 C2 )-1 (2-16)

where X3 is defined by

= 33 (t)= h(t)/ho (2-17)

in which h(t) and ho are the membrane thickness in the deformed and undeformed

configurations, respectively.

2.6.2 Solution of the Dynamic Equations

Using the discretized form of the principle of virtual work in the Lagrangian

description, the expression for the Green-Lagrange strain tensor, and the second Piola-

Kirchhoff stress tensor, the dynamic equation of motion for one triangular element in

global coordinates can be written as:

aI, = c5D. (MeDe + F~" Ft), (2-18)









where Ie is the virtual work for one element, 6De is the virtual nodal displacement

vector, Me is the mass matrix, De is the nodal acceleration vector, Fn"t is the internal

force vector due to membrane deformation, and FeX" is the external load vector.

Subscript e indicates values for one triangular element. The expressions for the element

mass matrix, Me, the internal force vector, F"t, and the external nodal load, Fet, can be

found in Lian (2003a). Equation (2-18) gives the virtual work for a single element. If we

sum all elements in global coordinates and take into account that the virtual work of the

whole system is zero, the governing equation for the membrane structure dynamic

response under external loads is:

MD(t) +Fnt = F- (2-19)

where M is the positive definite mass matrix (which remains constant) D(t) is the nodal

displacement vector for all nodes, D(t) is the nodal acceleration vector, Fint is the

internal elastic resisting force of the membrane due to its deformation, and F"" is the

nodal external load vector.

To solve the dynamic membrane equation (equation (2-19)), either an implicit or an

explicit scheme can be used for marching in time. The explicit scheme is computationally

economical for each time step and requires less storage. However, the time step needs to

be small in order to satisfy the stability requirements. The implicit scheme requires more

computational time and storage per time step because the system matrix is iteratively

updated at each time step to account for the nonlinear effects. But, generally, they are

unconditionally stable, and there is no restriction in the choice of time step size. In this

work, the implicit Wilson-O method (Wilson, 1973) is used. It is a one-step method which









is easier to code and behaves well in fluid-structure interaction problems. The Wilson-0

method has a dissipative character; therefore a smaller time step than that required by the

overall accuracy criterion is used to avoid excessive dissipation. A larger time step tends

to cause more numerical dumping as shown by Lian (2003a). Also, the method is

unconditionally stable for linear problems and is second-order accurate in time. For

details and illustrations of the membrane model and its computations results, see Lian et

al. (2002), and Lian (2003a).

2.7 Coupling between Fluid Solver and Structural Solver

During flight, the surrounding flow exerts aerodynamic forces on the membrane

structure which leads to deformations in the wing shape. At the same time, the change in

shape of the solid surface modifies the flow around the wing. To capture the interaction

between fluid and structure, the governing equations for both fluid and structure need to

be solved at the same time. However, structure and fluid formulations have different

mathematical and numerical properties, and depending on the problem complexity,

different techniques need to be developed to provide coupling between the computational

fluid dynamics (CFD) and computational structure dynamic (CSD) methods. One

approach is to treat the flow equations and the structural equations as one system of

equations (Bendiksen, 1991) and solve the system of equations simultaneously on a

single grid. However, the presence of multiple time scales in the fluid and structure

system makes a unified treatment inefficient in handling the computational stiffness.

To adequately handle different characteristics of the flow and structure interaction,

many researchers choose to solve the fluid and structure equations separately and couple

them through a synchronization process that involves subiterations (Smith and Shyy,

1997; Gordnier and Fithen, 2003; Zwaam and Pranata, 2002; Lian and Shyy, 2005; Lian









et al., 2002). The fluid solver and the structural solver run independently on their own

computational grid with their own choice of time step. Then a subiteration between the

two solvers is employed at each time step to couple the structure and fluid system of

equations and avoid phase lag.

Under the moving boundary framework for fluid-structure interaction problems, the

CFD grid needs to be automatically re-generated after the shape change. Since different

grids are used by the structure and fluid solvers, methods to transfer the information

between the CFD and CSD grids need to be employed (Smith and Hodges, 2000). Lian et

al. (2002), and Lian and Shyy (2005), in their study of the membrane wing dynamics,

have adopted a thin plate spline (TPS) interpolation method to exchange information

between the fluid and structure solver. The same method is also used in this work. The

thin plate spline interpolation is a global interpolation method which is invariant with

respect to rotation and translation, and thus is a good choice for interpolation on moving

surfaces. The interpolation function for the one-dimensional case can be written as:

N
H(x) ax x x -12logx- x, (2-20)
l=1

where H(x) is the interpolated displacement, a, is the undetermined coefficient, N is the

number of structural nodes on the interface, and x, is their location. Once the coefficients

a, are found based on the structural grid locations, the nodal displacement vector defined

on the CSD grid, D, is mapped to the CFD grid through an interpolation matrix, G,

X = GD, (2-21)

where 6X is the corresponding displacement on the CFD grid. Similarly, the pressure

forces can be mapped from the CFD grid to the CSD grid.










Read the CFD and CSD grids


Start the time loop, t = 0


Fluid Solver

Evaluate pressure p, by solving the Navier-
Stokes equations on the CFD grid.


Interface

Transfer the new pressure p, to the CSD grid
through TPS interpolation


4 Structural solver

2 Based on interpolated p evaluate the nodal
Displacement vector D,.


Interface

o Map the displacements from CSD to CFD
o grid using TPS interpolation: 6X,=GD,


Re-meshing

Based on the interpolated displacements
update the CFD grid


Re-compute the CFD grid metric terms

Update the Jacobian to satisfy the GCL


Next time step, t = t + At


STOP


Figure 2-5. Generic algorithm for fluid-structure interaction

In order to conduct time-accurate computations for fluid-membrane dynamics,

iterations need to be done at each time instant to enforce the fluid and structure coupling.

A generic algorithm for fluid-structure interaction problems with separate fluid and






50


structure solvers that exchange information between the CFD and CSD grid through

interpolation, and uses subiterations to enforce the fluid-structure coupling, is shown in

Fig. 2-5.














CHAPTER 3
FIXED WING MICRO AIR VEHICLES

3.1 Overview

Micro air vehicles with a maximum dimension of 15 cm or less and a flight speed

of 10-20 m/s are of interest to both military and civilian applications. Equipped with a

video camera or a sensor, these vehicles can perform surveillance and reconnaissance,

targeting, and bio-chemical sensing at a remote or otherwise hazardous location (Shyy et

al., 1999). With the rapid progress made in structural and material technologies, and

miniaturization of power plants, communication, visualization, and control devices,

several groups have developed successful MAVs (Pomsin-Sirirak et al., 2000; Grasmeyer

and Keennon, 2001; Ifju et al., 2002; Jones and Platzer, 2003; Jones et al., 2004).

With its reduced size, the MAV gains certain favorable scaling characteristics,

including low inertia and reduced stalling speed (defined as the minimum speed at which

the wing can produce sufficient lift force for flight) (Shyy et al., 1999). However, in

terms of aerodynamics, control, range and maneuverability, many outstanding issues

exist. MAVs fly in a low Reynolds number regime (104-105) due to their low flight speed

and small dimensions. This is often accompanied by laminar boundary layer separation,

transition, and reattachment, which leads to a decrease in the aerodynamic performance

(Shyy et al., 1999). Valuable information in these regards is reviewed by Lissaman

(1983) and Tani (1964), as well as addressed in several proceedings by Mueller (1998,

2000). Active control approaches with micro-electronic-mechanical systems (MEMS) for

MAVs are discussed by Gad-el-Hak (2001).









In the development of MAVs, there are three main approaches: (1) using flapping

wings for generating lift and thrust; (2) employing rotating wings; (3) using fixed wings

for generating lift. From the design point of view, fixed wing configuration is simple in

concept and easy to implement. On the other hand, flapping wings, which show

superiority over fixed wings in terms of lift generation and maneuverability at low

Reynolds number, is unfortunately difficult to put into practice due to complex wing

kinematics and power requirements. The simplicity in design of the fixed wing MAV, as

a miniature of large airplane wings, is counterbalanced by a deterioration in performance

because of the presence of separation bubbles, as its operating Reynolds number drops to

the range of 104 to 105. Also, the efficiency of conventional propellers diminishes with

decreasing diameter. Flapping wing design, which is predominantly selected as the

optimal approach by nature, compensates the complex kinematics by simultaneously

providing both lift and propulsion.

Figure 3-1(A) demonstrates a fixed wing design by Ifju and co-workers at the

University of Florida (Ifju et al., 2002). This design features a unique flexible wing,

which is advantageous to the MAV flight because it can adapt to wind gusts, offering a

smoother platform. Also, experimental studies show that the flexible wing has higher stall

angles than its counterpart rigid wing. The membrane wing consists of a leading edge

spar and six chord-wise battens, which are made of unidirectional carbon fiber prepreg

laminates. The membrane material is bonded to the spar and battens. While the

membrane is flexible, and cannot sustain bending moments, the carbon fiber is rigid, can

sustain bending moments and provides structural support.









The other approach, inspired by nature, uses flapping wings to produce both lift and

thrust. This concept has been practiced in the designs of Pornsin-Sirirak et al. (2000) and

is shown in Figure 3-1(B). A mixed design is presented by Jones et al. (2004). Their

MAV has a fixed wing used to generate the lift force, while using small flapping wings to

generate thrust.


Figure 3-1. Fixed and flapping wing MAVs. A) University of Florida 15 cm (6") MAV
design with a flexible wing. B) CalTech/AeroVironment Microbat flapping
wing MAV.

In this chapter we focus on the aerodynamics of both rigid and flexible fixed wing

MAVs, based on the designs of Ifju's group at the University of Florida. The steady, rigid

wing aerodynamics is presented first, followed by the unsteady flow over a flexible wing.









3.2 Rigid Wing Micro Air Vehicles

3.2.1 Wing Geometry

Figure 3-2 shows the wing shapes designed by Ifju and co-workers for their MAVs

and are used as a base for the present numerical simulations. Figure 3-2(A) shows an

older design of a single curvature airfoil wing, which has a 15 cm (6") span, a root chord

of 13.3 cm (5.24") and an area of 160 cm2 (24.8in2). A more recent design is shown in

Figure 3-2(B). The wing has a span of 12.7 cm (5"), a root chord of 10.4 cm (4.09") and

an area of 110.4 cm2 (17.11in2). The wing has two distinct regions: a central panel with

almost constant airfoil profile and no dihedral angle; and an outer panel with variable

cross-section and a dihedral angle of almost 7. The double curvature airfoil cross-section

offers better transversal (pitch) stability characteristics for flying wing airfoils.

A) I B)Z












Figure 3-2. MAV's wing shapes. A) 15cm (6") MAV wing shape. B) 12.7 cm (5") MAV
wing shape.

3.2.2 Computational Setup

To compute the aerodynamic performance of the MAV wing, the steady state,

Navier-Stokes equations for the incompressible flow are solved as described in Chapter

2. A pressure-based algorithm is used to solve the three-dimensional Navier-Stokes

equations written in curvilinear coordinates. Both first-order and second order upwind









schemes are employed for convection terms and second-order central difference schemes

are adopted for pressure and viscous terms.

For the flow over the rigid MAV wings, the physical problem under consideration

is the steady state flow over a wing in an unbounded domain. For all simulations, a

structured, multi-block grid is employed in the computations. The outer boundaries are

placed far away from the wing to reduce the boundary influence on the flow near the

wing. The schematic of the computational domains for the 6" wing and 5" wing is

presented in Figure 3-3.

(-6c, 5c, -6c) J 4 seC)
(-6 5-6 OUTLET -.",, ~~, ^ OUTLET .







c -5c, 6c) (c


A INLET (9c, -5c, 6c) B INLET (1., .,4e)


Figure 3-3. Computational domain for the MAV wing. A) 15 cm (6") wing. B) 12.7 cm
(5") wing.

In this work and in similar studies (Lian and Shyy, 2005; Lian et al., 2002; Lian et

al., 2003c), only the wing is taken into consideration, and no effects from the fuselage

and propeller are accounted for. Since the freestream velocity is along the chord-wise

direction and no propeller is modeled, only half of the wing is modeled (based on the

symmetry condition). A uniform velocity is specified at the inlet boundaries, while for

the outlet boundaries a zero-gradient condition is imposed. On the wing surface, the

noslip condition is specified.









3.2.3 Tip Vortex Effect on the Rigid Wing Aerodynamics

On a finite wing, the pressure difference between the high-pressure region below

the wing surface and the low-pressure region above the wing generates a circulatory

motion around the wing tip that affects the wing aerodynamics. The total drag coefficient

for a finite wing at subsonic speed can be written as (Anderson, 1989):


CD = CD,P +CD,F + L
T-e -AR (3-1)

where CD,P is the drag coefficient due to pressure, CD,F is the drag coefficient due to skin

friction, e is the span efficiency factor (which is less than one), AR is the aspect ratio, and

c L
C D 1 -- --
D, -e-.AR is the induced drag coefficient due to the existence of tip vortices.

Equation (3.1) demonstrates that the induced drag varies as the square of the lift

coefficient; at high angles of attack the induced drag can constitute a substantial portion

of the total drag. The MAV wings presented by the UF MAV group (Ifju et al., 2001; Ifju

et al., 2002; Albertani et al., 2004), seen in Figure 3-2, have a low aspect ratio of

approximately 1.4; the wing shape is dictated by the need to maximize the wing area, and

hence the lift, for a given dimension. Equation (3.1) illustrates that as AR is decreased,

the induced drag increases. Therefore, for these low-aspect wings, it is important to

investigate the tip vortex effects on the wing aerodynamics. In general, tip vortex effects

are two-folded:

* Tip vortex causes downwash that decreases the effective angle of attack and
increases the drag force (Anderson, 1989).

* Tip vortex forms a low-pressure region on the top surface of the wing, which
provides additional lift force (Mueller and DeLaurier, 2003).









Lian et al. (2003c), and Lian and Shyy (2005) performed numerous computations

to study the incompressible, laminar flow over the rigid 15 cm (6") MAV wing (Fig. 3-

2(A)). The relative position of the wing with respect to the outer-boundaries of the

computational domain is presented in Figure 3-3(A). Based on the root chord length and

a free-stream uniform velocity of 10m/s, the Reynolds number is 9x104. They studied the

flow over the wing at angles of attack of 6, 150, 270, 39, and 510. Figure 3-4 shows the

vorticity contours and pressure on the wing surface as they evolve from 6 to 270 angle of

attack. At 6 angle of attack, the tip vortices are clearly visible, but their strength is

modest and affects only a small portion of the wing. As the angle of attack increases to

270, the tip vortices are stronger and the low-pressure area on the wing's upper surface

becomes larger.
















A) B)

Figure 3-4. Streamlines, vorticity and pressure on the 6" MAV wing's upper surface as a
function of the angle of attack. A) 6 angle of attack. B) 270 angle of attack.
Adopted from Lian et al. (2003c).

At higher angle of attack (450), the flow separates from the leading edge, and the

low-pressure region associated with the tip vortex core increases, which helps to maintain

the increase in lift force.









Lian et al. (2003c) and Lian and Shyy (2005) concluded that for low-aspect ratio

wings, the tip vortex has an important effect on the wing aerodynamics. On one hand, the

tip vortex increases the lift at high angle of attack, similar to the vortex lift phenomenon

observed on delta wings (Torres and Mueller, 2001). On the other hand, the tip vortex

increases the induced drag, and therefore decreases the aerodynamic performance.

If we look at the MAV as a surveillance platform, it is desirable to have good flight

efficiency, which translates to larger values of the lift-to-drag ratio. To accomplish this

objective we try to maintain/increase the lift while decreasing the induced-drag. Various

methods to reduce the induced drag by decreasing the tip vortex effects are described in

the literature and confirmed by actual applications to aircraft wing design (LaRoche and

Paffy, 1996). In this study, the focus is on the implication of placing endplates at the

wing tip, as this is the simplest method from a manufacturing point of view.

In this work, we study the wing tip vortex effect on the rigid 15 cm (6") MAV wing

(Fig. 3.2(A)) and the influence of an endplate on the overall flow. Since a curved

endplate deteriorates the overall aerodynamic performances (Viieru et al., 2003), a

modified wing platform that allows straight endplates was created. Three wing

geometries have been studied: (1) the original wing (Fig. 3-5(A)) with dimensions

specified in section 3.2.1; (2) the modified wing (Fig. 3-5(B)) with a span of 14 cm (5.5")

and a wing area of 155 cm2 (24 in2), while the root chord has the same length as the

original wing; (3) the modified wing with endplates attached. The endplate has a length

of 4.4 cm (1.73") and a height of 3.4 cm (1.34"). The modified wing with endplates is

shown in Fig. 3-5(C).









Original wing shape Modified MAV wing Modified MAV wing and Endplates






A B C

Figure 3-5. Wing shape geometry. A) Baseline wing. B) Modified wing. C) Endplates
location on the modified wing. The arrow indicates the flow direction.

The computational domain is identical to that presented in Fig. 3-3(A). The

boundary conditions are described in section 3.2.2, with the addition of noslip boundary

conditions assigned to the endplate surfaces. The free-stream velocity is 10 m/s for all

cases, which gives a Reynolds number of 9x104 based on the root chord length. Lian and

Shyy (2005) have performed computations for the same original wing geometry (Fig. 3-

5(A)) using three different grid sizes ranging from 180,000 nodes (41x31 on the wing

surface) at the coarse level to 2.3 million nodes (101x61 on the wing surface) at the fine

level. They concluded that the coarse grid solution under-predicts the lift-to-drag ratio

(CL/CD) by 0.5% and form drag force by less than 0.7%, compared to the finest grid

solution. A similar result was also obtained for the lift force. However the skin friction

drag is sensitive to grid density, and the computation based on the coarse mesh

overestimates it by 10%. Since in the presently considered parameter range, the form

drag is at least one order of magnitude higher than the viscous drag, we use the coarse

grid in our numerical simulation to save computational time.

In Fig. 3-6, the vorticity magnitude contours behind the wing are plotted in planes

perpendicular to the x-axis, i.e., the stream-wise direction. The streamlines are also

displayed to clearly show the wing tip vortex. It can be seen that both wings without the









endplate show a strong wing tip vortex, while for the modified wing with endplates, the

tip vortex core is smaller and dissipates faster.


Figure 3-6. Vorticity magnitude contours behind the wing at 6 angle of attack. A)
Original wing. B) Modified wing. C) Modified wing with endplates. Arrow
indicates the flow direction.

The endplate basically deters the flow from the high-pressure region on the lower

wing surface from reaching the low-pressure region on the upper wing surface. In the

absence of the endplate, the fluid is accelerated towards the wing tip by the vortex core's

low-pressure zone present on the upper wing surface, as indicated by the streamlines

plotted in Fig. 3-7(A), (B). The low-pressure zone associated with the tip vortex core can

be clearly identified in the plot of the span-wise pressure variation on the wing's upper

surface (Fig. 3-9(A) and Fig. 3-10(A)).










Wing tip

--

---





: p -90-79-68-57-46-35-24-13 -2 9 20


B
Wing tip






as_


r


Wng tip

--._ --


K


p: -90 -79-68 -57-46 -35-24 -13 -2 9 20


Figure 3-7. Pressure contours on the wing's upper surface and streamlines slightly above
the wing at 6 angle of attack. A) Original wing. B) Modified wing. C)
Modified wing with endplates. Pressure is in [N/m2].

The bending of the streamlines towards the wing tip in the absence of the endplates,

caused by the pressure difference between the high-pressure zone (below the wing

surface) and the low-pressure zone (above the wing surface) can also be observed by

looking at the streamlines slightly below the wing surface, as shown in Fig. 3-8(A) and

(B). One may also note that the tip vortex is stronger in the case of the modified wing

without endplates than in the case of the original wing. This points to a more significant

pressure drop on the lower wing surface, and can be observed in Fig. 3-8(B). The

decrease of the high-pressure area due to the circulatory motion in the absence of the


p -90 -79 -68-57-46 -35 -24 -13 -2 9 20


I


I _







62


endplate can also be recognized by looking at the span-wise pressure coefficient on the

lower wing surface, which is plotted in Fig. 3-9(B) and Fig. 3-10(B).


|__ VWng tip

A p: 30-24-18-12-6 0 6 12 18 24 30
p -30-24-18-12 -6 0 6 12 18 24 30


LW Vng tip


p -30-24-18-12-6 0 6 12 18 24 30


SVWng tip


p: -30 -24 -18 -12 -6 0 6 12 18 24 30


Figure 3-8. Pressure contours on the wing's lower surface and streamlines slightly below
the wing at 6 angle of attack. A) Original wing. B) Modified wing. C)
Modified wing with endplates. Pressure is in [N/m2].

The endplate decelerates the flow near the wing tip. On one hand, this decrease in

velocity reduces the pressure drop on the upper wing surface corresponding to the vortex

core (Fig. 3-9(A) and Fig. 3-9(B)), reducing the vortex lift. On the other hand, a lower

velocity slightly below the wing increases the high-pressure area since more momentum

is transferred to the wing as pressure, instead of being shed as vorticity at the wing tip.

The increase in the high-pressure zone on the lower wing surface in the presence of the

endplate can be clearly seen from the span-wise pressure coefficient on the lower wing








63



surface plot (Fig. 3-9(B) and Fig. 3-10(B)), and also from the pressure contours on the


lower wing surface shown in Fig. 3-8.


02-
0A)

A)


01 02 03
Z / Local span


0a
C-


04 05


Z / Local span


Figure 3-9. Pressure coefficient on the wing surface at x/c = 0.34 and 6 angle of attack.

A) Upper wing surface. B) Lower wing surface.


-e- Original wing, NO endplates
-v-- Modified wing, NO endplates
--,-- Modified wing, WITH endplates


0 15 4 I
-e- Original wing, NO endplates
01 --v-- Modified wing, NO endplates
S05 --0-- Modified wing, WITH endplates
0 ------------- ----------
0--
9-0 05 --
-01---
-015- -
-02------ ------


01 02 03
Z / Local span


1 -0 :
04 05


01 02 03
Z/ Local span


Figure 3-10. Pressure coefficient on the wing surface at x/c = 0.53 and 6 angle of attack.
A) Upper wing surface. B) Lower wing surface.


To get a better picture of the endplate effects, the pressure difference between


lower and upper wing surfaces was integrated along the local chord direction at different


span-wise locations. The result is the lift and drag estimation in the span-wise


distribution. Figure 3-11 shows the span-wise lift distribution. It can be noticed that the


endplate, by blocking the flow around the tip, reduces the downwash. Therefore the


--e- Original wing, NO endplates
S- -v-- Modified wing, NO endplates
--A-- Modified wing, WITH endplates
-0 1 0

-02 ---- -






-0 -

- 0 o01 012 013 04 05


02A
0A)

A)


04 05










effective angle of attack is increased, leading to an increase in lift towards the wing tip.

The figure shows that when the endplates are attached, the lift is increased when

compared with the wing shape without the endplate.

45




35
4 5-------------,------------------





25 -

-e- Original wing, NO endplates
2 *- Modified wing, NO endplates
S- Modified wing, WITH endplates

0 005 01 015 02 025 03 035 04 045 05
Z / Original Wing Span

Figure 3-11. Span-wise lift distribution at 6 angle of attack.

The modified wing with endplates produces almost the same amount of lift as the

original wing. However, if we look at the span-wise drag distribution plotted in Fig. 3-12,

the modified wing with or without endplate offers the lowest drag for almost 75% of the

half wingspan (starting from the root).

After computing the total lift and drag by integrating the pressure difference on the

wing surface, we conclude that for an angle of attack of 6, an improvement in lift-to-

drag ratio of 10.1% for the modified wing configuration with endplates is observed as

compared to the baseline configuration (original wing, no endplate) (Viieru et al., 2004).

This improvement is based mainly on the reduction in drag offered by the modified wing

shape, since the lift force profiles are similar.










For 150 angle of attack, the modified wing with endplates shows an increase of

1.4% in lift-to-drag ratio compared to the baseline configuration. This is achieved by an

important reduction of drag, while lift is decreasing with by 5%.



-e- Original wing, NO endplates
055 -v Modified wing, NO endplates
-A Modified wing, WITH endplates


045 .













3.3 Flexible Wing Micro Air Vehicles
035- 'a






As mentioned before, flexible wings increase the aerodynamic performance of
025
0 005 01 015 02 025 03 035 04 045 05
Z /Original Wing Span

Figure 3-12. Span-wise drag distribution at 60 angle of attack.

3.3 Flexible Wing Micro Air Vehicles

As mentioned before, flexible wings increase the aerodynamic performance of

fixed and flapping wings MAVs, but due to the complexities involved in the fluid-

structure problem, less rigorous work has been considered in this regard. Waszak et al.

(2001) performed an experimental study of a flexible wing MAV. Lian et al. (2003c)

numerically simulated the flow over a flexible wing. They showed that a MAV with a

flexible wing yields better performance than that of a rigid wing because of its shape

adaptation and delay to stall. Ho et al. (2003) experimentally demonstrated that flexible

wings have better lift production for flapping flight. However, no decisive conclusion has

been drawn on how the wing flexibility enhances the aerodynamic performance.

The influence of the wing flexibility of a 12.7 cm (5") MAV at 60 angle of attack is

studied next. The geometry is identical with that presented in section 3.2.1 (Fig. 3-2(B)).










The central wing panel is rigid, while the outer panels are covered with a flexible material

that has a uniform thickness of 0.2 mm and a density of 1200 kg/m3 (Ifju et al., 2001; Ifju

et al., 2002; Lian et al., 2002). The flexible part has a rigid leading edge spar and two

flexible carbon-fiber battens, as shown in Fig. 3-13(B). The battens are modeled (in the

structural solver) as membrane material with increased density to simulate the dynamic

response. For the finite element structural solver (chapter 2), an unstructured mesh with

598 nodes and 1076 elements was created. A coarse CFD grid with 210,000 nodes was

used, since the overall aerodynamic performances are not very sensitive to grid size (Lian

and Shyy, 2005) and significant computational time can be saved. Both CFD and FEM

grids are presented in Fig. 3-13.













Grey area = Rgid part,
Red = Wing tip regions; Green = Rigid leading edge boundary;
Black = Rigid central panel Blue = Batten locations

Figure 3-13. Surface grid for the 5" MAV wing. A) CFD coarse mesh. B) FEM mesh for
the flexible part of the wing indicating the boundary conditions for the
structural solver.

Based on a free-stream uniform velocity of 10 m/s and the root chord length, the

Reynolds number for this case is 7.1x104. For the unsteady computations, a time step of

2x10-4 s is used to ensure that membrane vibration time scales are captured. Two sub-










iterations were used to ensure the time synchronization as shown by Lian et al. (2003c).

The fluid-structure framework is presented in chapter 2.

Even though the flow stream is uniform, the membrane exhibits self-initiated

vibrations, as one can see in Fig. 3-14. This was also observed by Lian et al. (2003c) in

their computations of the flow over a 15 cm (6") MAV with flexible wing.

Vertical dispacement 5in MAV
02


0 15 -- -------- -1 --r- ----I------ -


I0 15-- -- -- --- -I- - -
015




01





-o005-v-------- -- near wingtip
--- between the battens
-oil 1 L---
0 001 002 003 004 005 006 007 008 009 01
time [sec]


Figure 3-14. History of the vertical displacements of two points, one near the wing tip
and the other one between the two battens.

It can be clearly seen that the point near the tip exhibits larger displacements than

the point between the two battens. The probable explanation of this fact is the presence of

the strong circulatory motion of the tip vortex. Also, the battens' rigidity restrains the

movement of the points between them.

The vertical displacement contours are visualized in Fig. 3-15 at different time

instants. The first snapshot shows how the membrane starts to inflate in the region of

largest pressure difference. The ripples in the membrane, especially between the battens,

can be observed in the next time snapshot. This is correlated with the periodic behavior

seen in Fig. 3-14.











tme =0005 *,. time=015 00 le025 a

J1 U:::^ Ju^
4 ......
&L.. : & L


-i .B0045


., 0056
Y


u ?


Figure 3-15. Vertical displacement contours at various times. The displacement is in
[cm].


To assess the overall aerodynamic performances, the pressure and shear stress are


integrated over the wing surface at each time step. In Fig. 3-16 the lift and drag


coefficients history are presented for the 12.7 cm (5") MAV wing at 6 angle of attack.


Lift coefficient history, 5in MAV, AoA= 6 deg Drag coefficient history, 5in MAV, AoA= 6 deg


0075

S 007

I ffl 0 065

8 006

0 055

S 005


I I I


015 -- I I I -- --- 0045 --- I --1-----

0 001 002 003 004 005 006 007 008 009 0 0 001 002 003 004 005 006 007 008 009 01
time [sec] time [sec]

A) B)


Figure 3-16. Aerodynamic forces history for the 5" flexible wing MAV. A) Lift
coefficient. B) Drag coefficient.


l.= 0036


'4


;I,
nnz,
_u









Next, the time averaged aerodynamic performance is compared with the steady-

state, rigid wing at the same angle of attack and with experimental values measured in the

wind tunnel for the same rigid wing geometry (Viieru et al., 2004). The results are

compared in Table 3-1.

Table 3-1. Aerodynamic performances for flexible and rigid wing.
Parameter Computational Computational Experimental -
Flexible wing Rigid wing Rigid wing
(time average)
Lift coefficient (CL) 0.37 0.40 0.42
Drag coefficient (CD) 0.064 0.076 0.085
Lift / Drag ratio (CL/CD) 5.78 5.26 5.06


In the computation of the lift and drag coefficients, both pressure and viscous

forces are taken into consideration. The flexible wing's averaged lift coefficient is

smaller, but comparable with the rigid wing's computational and experimental values.

This indicates that for small angles of attack, the aerodynamic performances are not very

different. In terms of drag coefficient, it was found to be smaller for the flexible wing,

leading to a better lift-to-drag ratio. The best explanation for the drag reduction for the

flexible wing is the adaptation to the flow, resulting most probably in a reduction of the

tip vortex effects and therefore a reduction in the induced drag. The small difference

between the aerodynamic performance of the rigid and flexible wing was also noticed by

Lian et al. (2003c) in his computations for the 6" MAV wing. However, experimental

work has proven that the advantages of flexible wings come into play at high angles of

attack. Generally, membrane wings show a higher stall angle than rigid wings (Waszak et

al., 2001; Albertani et al., 2004). Due to complicated flow behavior at high angles of

attack for this range of Reynolds numbers, mainly flow separation and transition from









laminar to turbulent boundary layer, numerical investigation of the behavior of the

flexible wings at high angles of attack remains an open topic.

3.4 Summary and Conclusions for Fixed Wing Aerodynamics

In this chapter, we numerically investigated the flow over rigid and flexible fixed

wings. The steady-state, laminar, incompressible Navier-Stokes equations were solved

using a pressure based solver to study the flow over 6" and 5" wingspan wings at

Reynolds numbers ranging from 7x104 to 9x104 at different angles of attack. For the

flexible wing computations, the unsteady Navier-Stokes solver was coupled with a finite

element structural solver to investigate the fluid-structure interaction problem of a

flexible wing in a uniform flow.

First, the effect of the tip vortex on the overall aerodynamics of a low aspect ratio

15 cm (6") MAV wing at 60 and 15 angles of attack was studied. To reduce the induced

drag generated by the tip vortex, the endplate effect was also studied for both the same

wing configuration and for a modified wing that allows for the installation of a straight

endplate.

The following observations were made:

* The modified wing shape reduces the drag, with or without endplate. However the
lift is also reduced due to a stronger vortex and a reduced area.

* The endplate parallel to the freestream velocity increases lift without a significant
increase in drag. The increase in lift is closer to the tip, where the endplate reduces
the flow travelling from the lower surface to the upper surface of the wing,
reducing the downwash, and increasing the effective angle of attack. Consequently
the lift increases.

* The effectiveness of the endplate diminishes as the angle of attack increases, due to
a stronger wing tip vortex.

Secondly, the impact of wing flexibility was studied on a 12.7 cm (5") wingspan

wing at 6 angle of attack and Reynolds number of 7.1x104. The results showed that the






71


membrane wing vibrates, even in a steady free-stream. For small to moderate angles of

attack, and before stall limit, the time-averaged lift and drag coefficients of the membrane

wing seem to be comparable with those of a rigid wing. However, the flow behavior at

high angles of attack for this range of Reynolds numbers, which departs form the pure

laminar assumption, and from the high Reynolds number turbulent models, presents

increased difficulties for the numerical investigation of the flexible wings and remains an

open topic for research.














CHAPTER 4
AERODYNAMICS OF LOW REYNOLDS NUMBER HOVERING AIRFOILS

Even though three-dimensional effects are critical for MAVs, experiments and

computations in two-dimensions continue to provide important insight into flapping wing

aerodynamics. From the computational point of view, 2-D simulations require less

computational power, and grid sensitivity and convergence studies can be performed

more readily (compared to 3-D cases). Wang (2000a, 2000b) demonstrated that the

aerodynamics of 2-D wings could help explain some of the enhanced lift observed in

insect flight. In this chapter, the flow structures around a hovering 2-D elliptic airfoil are

studied. First, the viscous force computation is validated for the Reynolds number range

encountered in insect flight. Second, the results for a hovering airfoil are compared with

similar numerical and experimental work to validate the current fluid solver. Finally, the

lift generation mechanisms are compared for two hovering modes, and the effect of the

Reynolds number on the aerodynamics is investigated.

4.1 Validation of Viscous Force Computation

Low Reynolds number flows are characterized by large viscous forces that become

closer in magnitude to the pressure forces, and therefore have a large contribution to the

overall lift and drag forces. Generally, the shear stress at the wall is computed as the

product of the velocity gradient in the normal direction at the wall surface and the

viscosity:


Z- wall = -l (4-1)
On) all









In the present approach, equation 4-1 is discretized using a first order scheme:


Twal = (Utan )p (tan)wall (4-2)
An

where, (Utan)P is the tangential velocity at the center of the first cell from the wall,


(Utan)waii is the tangential velocity of the wall (non-zero for moving boundary problems)

and An is the normal distance from the wall to the center of the first cell near the wing

surface.

To validate the present formulation, two cases are studied. First, the steady flow

over a flat plate at zero angle of attack at different Reynolds numbers representative of

insect flapping flight is computed. The evaluated friction coefficient is compared with

analytical results. Secondly, the flow over an oscillating flat plate in a quiescent medium

is simulated in order to assess the viscous force computation for moving walls.

4.1.1 Steady Flow over a Flat Plate at Zero Angle of Attack

To simulate the steady flow over a thin plate of chord c = 1.0, the Navier-Stokes

equations are solved on a grid with 100 points along the streamwise and 60 points along

the cross-stream direction. The distance from the wall to the first cell center is 5x10-4c,

which guarantees a sufficient number of points in the boundary layer for high Reynolds

number flows. The numerical simulations are performed for Reynolds numbers (based on

chord length and free-stream velocity), of 10, 102, 103, and 104. The second order upwind

scheme was employed for convection terms, while second order central difference

schemes are adopted for the pressure and viscous terms. A no-slip boundary condition is

imposed on the airfoil surface.









In Fig. 4-1, the numerical and analytical velocity distribution in the boundary layer

for a flat plate at different Reynolds numbers is plotted. The analytical velocity

distribution is defined as (Schlichting, 1979):


u = Uf'(r)


(4-3)


where Um is the free-stream velocity, f'(q) is the solution of Blasius's ordinary


differential

coordinate.


equation (Blasius, 1908), and 7 = yVU /(v. x) is a dimensionless


U velocity profile, Re = 10,000


U velocity profile, Re = 100


Figure 4-1. Numerical and analytical velocity distribution in the boundary layer along
different locations on a flat plate and different Reynolds numbers. A) Re=104,
B) Re=103, C) Re=102. Here r = y Um /(v x) is the dimensionless vertical
coordinate.









Figure 4-1 shows good agreement between the numerical results and the analytical

formulation, especially for high Reynolds numbers. As the Reynolds number decreases,

the computed velocity distribution moves further away from the analytical distribution.

This is expected, since the Blasius's equation is based on Prandtl's (1904) boundary layer

equations, which implied the high Reynolds number flow assumption. The analysis of the

flow over a flat plate using Blasius' equation is restricted to semi-infinite plate, because

of the parabolic nature of the Prandtl's boundary layer equations that cannot account for

upstream changes in shear stress initiated by the trailing-edge of a finite-length plate. The

discrepancy between the analytical and numerical velocity distribution is visible at

locations near the trailing-edge of the plate, as one can observe in Fig. 4-1 for all

Reynolds numbers. Figure 4-2 also shows that as the Reynolds number decreases, the

friction coefficient given by Blasius departs from the numerical solution. To take into

account second order effects, generalizations of Prandtl's boundary layer equations were

developed. For the case of the flat plate, Stewartson (1974) and Messiter (1970) found an

expression that improves the prediction of the skin-friction coefficient for low Reynolds

number. The analytical expressions for skin-friction coefficient, defined for 2-D cases as

Cf = Drag/(0.5pcU2) are presented in Table 4-1.

Table 4-1. Numerical and analytical friction coefficient values for flat plate.
Case Reynolds Cf Blasius (1904) Cf Messiter (1970) Cf present
number c, =1.328/J/R C = 1.328//Re +2.668/(Re,)78 numerical
results
1 104 0.0133 0.0141 0.0142
2 103 0.0420 0.0483 0.0486
3 102 0.133 0.180 0.182
4 10 0.420 0.776 0.785











Figure 4-2 shows a very good agreement between the improved analytical solution


of Messiter (1970) and the numerical results for a large range of Reynolds numbers,


therefore validating the viscous force computation method employed in the solver.


Friction coefficient vs Reynolds number


06--
..-.---- -------------- Blasius
0 4 1 1 1 1 1 1 I
04" -- -- Messiter

02 Stream


01 4-
S I- -IrTTTr P I I IT* I J T rr
I TTr I IT S rF I F T T- TT F
006
I--- I I I-4 I -
0 04 I I -


0 02
002 -2-- I I ++ 1 4 + I 4 -


001 L L L I I---------------
101 102 103 104
Re


Figure 4-2. Numerical and analytical friction coefficient Cffor a flat plate versus
Reynolds number.

4.1.2 Flow over an Oscillating Flat Plate in a Quiescent Environment

In the previous sub-section, the viscous force computation is validated for steady


flows at Reynolds numbers from 10 to 104. Next, the concept needs to be validated for


the case of moving walls (wall velocity is non-zero). If we consider an infinite plate


oscillating along the x-axis, the wall velocity is defined as:


u(y = 0, t) = Uo cos(ot) (4-4)


where Uo is the maximum velocity and c) = 27rf ; being the oscillation frequency. The


integration of equation 4-4 gives the plate displacement:


h(t) = h, sin()t); I i it h = Uo l/ (4-5)











The analytical velocity field for this motion is given by Stokes (1851):


u U/U = exp(-r) cos(ot r/) (4-6)


where r7 = yJ-) /2v is a non-dimensional cross-stream coordinate.


Following the definition given in equation 4-1, the analytical wall shear stress for

an oscillating plate is given by (White, 1991):


Trwa = U0o p sin(t- r/4) (4-7)


The unsteady flow over an oscillating plate of length 1 is solved on a grid with 50

points along the chord and 80 points in the vertical direction using the Navier-Stokes

solver. The distance from the wall to the first cell center is 103.


Wall shear stress for oscillating plate
0.15
r computed
wall/1O computed
0 T analytic sol
0.1( - - -l-/10 analytic sol
0.05






0.05
0.1




-0.1


-0.15 -I I
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
time/T,


Figure 4-3. Skin-friction and wall velocity for an oscillating plate. Oscillating frequency
f = i/2ir =1/I r, amplitude ha=0.5 and Reynolds number =103. The wall
velocity is scaled to an order of magnitude comparable to the wall shear stress.
Continuous lines = present computation, Symbols = analytical solution for
infinite plate.










For the case studied, o) = 2 results in a frequency fof 1/7 and a period T of 7t. To

ensure a sufficient number of points for one oscillation cycle, a time step size of 0.01 is

chosen. Based on the maximum velocity Uo and chord length, the Reynolds number is

103.

In Fig. 4-3, the numerically predicted wall shear stress on the finite plate is

compared with analytical values for an infinite plate given by equation 4-7. The figure

shows a good agreement between numerical and analytical values. The small phase shift

between numerical and analytical results can be explained by the effects of the leading

and trailing edges of a finite plate. In Fig. 4-3 the scaled wall velocity is also plotted to

contrast the lag between maximum shear and maximum velocity.


6 1 -- ---------- ----- ----- -----------------


5 -- -







2 ------------- ----^- -- ,-- T = .1 -

2 -- - !- t/T 0.21 ----
t/T = 0.29 0.21

t/T 0.45 t/T = 0.05



-1 -0.75 -0.5 -0.25 0 0.25 0.5 0.75 1
ulU,


Figure 4-4. Velocity distribution above an oscillating plate. Oscillating frequency
f = ) 22 = 1/ i, amplitude ha=0.5 and Reynolds number =103. Continuous
lines = present computational results. Symbols = analytical solution for
infinite plate.









The theoretical velocity distribution above the plate, given by equation 4-6, is

plotted, along with the computed velocity field, in Fig. 4-4 for different time instants

during one half-period. The skin friction and velocity profiles presented in Fig. 4-3 and

Fig. 4-4 show a good agreement between theoretical and computed results.

Good agreement between the analytical solution for steady flow over a flat plate, as

well as for an oscillating plate is attained, validating the viscous force computation

methodology employed in the fluid solver for stationary or moving walls.

4.2 Kinematics of 2-D Airfoil in Hovering Flight

In hovering flight, insects and small birds draw the still ambient fluid and move it

downwards in a succession of vortex dipoles (or rings in 3-D) that carry the fluid

momentum and generate the lift necessary to stay aloft. The hovering flight plays an

important role for insects in feeding. Different hovering patterns are visualized in nature.

Figure 4-5 shows the figure "8" pattern observed in hovering hummingbirds.











Figure 4-5. Hummingbird in hovering flight. A) Forward stroke, B) Backward stroke.
Picture taken by Wei Shyy.

In an extensive survey, Weis-Fogh (1973) noticed that most insects, such as fruit

flies, bees, and beetles, employ symmetric back-and-forth strokes near a horizontal stroke

plane, as seen in Fig. 4-6(A). Weis-Fogh referred to this style as "normal hovering".

Some other insects, like true hoverflies and dragonflies employ asymmetric strokes along

an inclined plane. In this case the lift force is primarily generated during the downstroke










(Wang, 2005). The dragonfly's inclined plane hovering schematic is presented in Fig. 4-

6(B). Since the "normal hovering" is employed by most of the insects, it receives most of

the focus in experimental and computational investigations (Dickinson et al., 1999;

Ellington et al., 1996; Sun and Tang, 2002a, 2002b; Wang et al., 2004a). In this study, we

focus on two hovering modes, the "water treading" and the "normal" hovering modes.

The kinematics of these hovering modes are presented in the next section.

Normal Hovering Dragonfly
SLilt
Drag *Drag


raDrag


A B

Figure 4-6. The schematics of two hovering styles. A) Hovering using a horizontal plane.
B) Hovering using an inclined stroke plane. Sketched from Wang (2005).

4.2.1 Kinematics of "Water Treading" Mode

The "water treading" hovering mode is one of the cases that Liu and Kawachi

(1998) used to validate their finite volume algorithm. The hovering mode studied is based

on the so-called "water treading" mode as defined by Freymuth in his experiments

(Freymuth, 1990) and consists of a translational component (plunging), given by:

h(t) = ha sin(2kt + (p), (4-8)

and a rotational component (pitching), given by:

a(t) = ao + a, sin(2kt). (4-9)

In equation 4-8, ha is the plunging amplitude and k is the reduced frequency of the

sinusoidal oscillation. In equation 4-9, ao is the initial angle of attack, aa is the pitching









amplitude, and p is the phase difference between the plunging and pitching motion. The

reduced frequency k is defined by:

k = 2irfc/2Ue, = c/2ha (4-10)

where c is the airfoil chord length,f is the oscillation frequency, and Ur = 2- rfh is the

reference velocity (equal to the maximum plunging velocity). It should be noticed that the

reduced frequency, by definition, varies with the inverse of the stroke amplitude. By

setting the reference velocity (Uref) to one, the relationship between the reduced

frequency, k, and the frequency, f is: f=k/7c, and the stroke period is T=7c/k. The

Reynolds number is defined as Re= Ue c/v, where vis the kinematic viscosity of the

fluid. A schematic of the airfoil movement is shown in Fig. 4-7.


Y
Forward stroke




-ha ha


Backward stroke




-ha ha


Figure 4-7. Schematics of the "water treading" hovering mode.

The specific parameters for the airfoil movement will be specified for each

individual case in this study that employs "water treading" hovering mode.









4.2.2 Kinematics of "Normal" Hovering Mode

In this work a "normal" hovering mode based on a sinusoidal motion (Wang et al.,

2004) is also studied. The equations describing the elliptic airfoil motion are:

h(t) = ha sin(2kt + (p) (translation) (4-11)

a(t) = ao a sin(2kt) (rotation) (4-12)

where ha is the stroke amplitude, ao is the initial angle of attack, aa is the pitching

amplitude, ( is the phase difference between the plunging and pitching motion, and k is

the reduced frequency as defined in equation 4-10. A schematic of the motion is shown in

Fig. 4-8.


Forward stroke





X


Backward stroke



X

-ha7 ha

Figure 4-8. Schematics of "normal" hovering mode.

Based on the above equations, the parameters describing the airfoil motion will be

specified for each case studied in this work.

4.3 "Water Treading" Hovering Mode

Liu and Kawachi (1998) placed the airfoil in a body-centered inertial frame of

reference that undergoes the hovering motion about a fixed axis at origin. The plunging







83


motion is implemented in the boundary conditions at the open boundary stencils. The O-

type grid used in their work has different sizes, varying from 71x71 to 145x99 grid

points, in order to investigate the effect of the grid resolution. The minimum grid spacing

adjacent to the solid wall is 2.4x10-3c.

In this work, an O-type grid with different resolutions was also employed. At the

coarse level, the grid has 90x70 nodes, while at the fine level the grid has 250x120 nodes.

The computational domain for the coarse grid is shown in Fig. 4-9. The minimum grid

spacing adjacent to the airfoil varies from 2.4x10-3c (coarsest grid) to 1.0x10-3 (finer

grid).

DgrO (90x70, dy=2.4e-03)







>- o 4



-10 -

-10 0 10
x


Figure 4-9. The coarse multi-block grid used in the current work. The grid has 90 points
around the airfoil and 70 points in the radial direction. The distance from the
airfoil to the outer boundary is 15 chords.

The unsteady, laminar, incompressible, Navier-Stokes equations are solved using

the pressure based algorithm described in chapter 2. To handle the large airfoil

displacements and rotations, a moving grid algorithm is coupled with the flow solver.









4.3.1 Grid and Time Increment Sensitivity Analysis

The spatial accuracy of the present algorithm is examined by employing three grid

size levels. A refinement of the time step size was also performed on the intermediate

grid (named Dgrl). The hovering mode is described by equations 4-8 and 4-9 with the

same kinematics parameters as in the "water treading" hovering mode presented earlier.

Table 4-2 provides a summary of the cases run, along with grid and time increment

details. A schematic of the movement is presented in Fig. 4-7.

Table 4-2. Grid and time increment details for the "water treading" hovering cases. Here
6y represents the distance from the wall to the first grid point and 6t is the
time step size. Re=1,700, ha=l.0c, ao~=0, a,=420, k=0.5, p = -90.
Case Grid Computational time (hrs) Number of time
Case Grid 5t
(single processor) steps
1. DgrO (90x70, 6y=2.4x10-'c) 0.01 7.5 (Intel Itanium 1.3GHz) 5600 (~ 8 periods)
2. Dgrl (190x80, y =1.0x10-'c) 0.01 31.4 (Intel Itanium 1.3GHz) 5600 (~ 8 periods)
3. Dgr2 (250x120, 6y =1.0x10'c) 0.01 64.3 (Intel Itanium 1.3GHz) 5600 (~ 8 periods)
4. Dgrl (190x80, y =1.0x10-'c) 0.02 13.3 (Dec Alpha 850MHz) 2800 (~ 8 periods)
5. Dgrl (190x80, 6y =1.0x10-c) 0.005 44.8 (Intel Itanium 1.3GHz) 11200 (-8 periods)


The force coefficients are obtained by normalizing the forces by: 0.5pV(t)2 where

p is the density and V(t) = 0.5U2 (Liu and Kawachi, 1998).

Figure 4-10 and 4-11 show the lift and drag coefficient history for three different

grid resolutions. The results show a good agreement between solutions and indicate that

reasonable grid independence was obtained. The lift and drag forces attain periodic

behavior after almost 5 periods. The highest lift peak values are observed for the finest

grid (250x120).

To study the effect of time step size with a good spatial accuracy, the intermediate

grid, named Dgrl (190x80), is used. Figures 4-12 and 4-13 show the lift and drag

coefficient history for time increments varying from 0.005 to 0.02











Lift coefficient history for St=0.01


....... DgrO ( 90x70)
--- Dgrl (190x80)
-Dgr2 (250x120)


time/T


Figure 4-10. Lift coefficient for three periods and different grid sizes and St=0.01.
ha=i.0c, ao=0o, aa=42, k=0.5 and Re=1,700.


Drag coeffcient history for St=0.01


6.5
time/T


....... DgrO ( 90x70)
--- Dgrl (190x80)
- Dgr2(250x120:
7.5 8


Figure 4-11. Drag coefficient for three periods and different grid sizes and tt=0.01.
ha=1.0c, ao=0o, aa=42, k=0.5 and Re=1,700.










Lift coeffcient history on grid Dgrl (190x80)


Figure 4-12. Lift coefficient for three periods and different time increments. ha=1.0c,
ao=0, aa=42, k=0.5 and Re= 1,700.


Drag coeffcient history on grid Dgrl (190x80)


6.5
time/T


Figure 4-13. Drag coefficient for three periods and different time increments. ha=l.0c,
ao=0, aa=42, k=0.5 and Re=1,700.









As observed in Fig. 4-10, the lift coefficient reaches a plateau close to the

maximum pitch angle (the middle of forward/backward stroke), and then a sharp increase

in lift is noticed as the airfoil starts pitching down. This is captured well for time

increments of 0.005 and 0.01; a smaller time step better captures the vortex shedding

phenomenon responsible for the sudden decrease in lift. Using a time step size of 0.02,

the sharp increase in lift following the pitch-down of the airfoil is not captured very well.

Based on the results for different spatial and temporal resolutions, we conclude that

the combination of an intermediate grid (190x80) and a time step size of 0.01 gives

satisfactory balance between numerical accuracy, the flow structures, and computational

cost. The results also indicate that reasonable grid independence was obtained, taking into

account the unsteadiness of the flapping flight.

4.3.2 Comparison with Similar Work

Based on the grid and time step size refinement study, the results obtained on the

intermediate grid (190x80) with a time increment of 0.01 are compared with the work of

Liu and Kawachi (1998). Figure 4-14(A) shows the lift and drag coefficients history for

one period as obtained by Liu and Kawachi, while Fig. 4-14(B) shows the current

computed values. The lift and drag coefficients obtained in the present numerical

simulation correspond well to those of Liu and Kawachi, and further confirm the

capabilities of the Navier-Stokes solver used in this work.

4.3.3 Lift Generation Mechanism for "Water Treading" Hovering Mode.

In this section, the lift generation mechanism is investigated for this particular

hovering mode. The flow analysis is based on the last computed period in order to avoid

the transient effects and to allow the force variation to reach a periodic pattern. The data










obtained using the finest grid (250x120) and a time increment of 1x10-2 is used for this

analysis. Eight time instants are selected for analysis during one hovering period and are

identified in Fig. 4-15. The airfoil position and pitching angle, as well as the translational

and rotational velocities, are plotted in Fig. 4-15(A). Fig. 4-15(B) shows the lift and drag

coefficient variations.

Resolution of time increment Lift and Drag coefficients for one period



S3 ... ...... C




I L
-42



-4 --19 20 21 22 23 24 25
19 20 21 22 23 24 25 Time
Time

Figure 4-14. Lift and Drag coefficients for one period. h,=1.0c, ao=0, a,=42, k=0.5 and
Re=1,700. A) Liu and Kawachi (1998) (figure obtained via personal
communication). B) Present computational results.

At the beginning of the forward stroke, the airfoil accelerates from rest and starts to

pitch-up rapidly. Fig.4-15(B) shows that most of the lift and drag force is generated when

the airfoil travels from time instant A to B. Between time instant B and C, the airfoil

travels close to the maximum translation velocity and at an almost constant maximum

pitching angle (Fig. 4-15(A)). Figure 4-15(B) shows that, during this period, the lift curve

reaches a plateau.

The clockwise vortex shed earlier (Fig. 4-16(A) at time A) accelerates the flow

over the upper surface of the airfoil, and creates a low-pressure zone (Fig. 4-16(A) and

(B), time B).








89



Kinematics for "water treading" hovering mode
12
1 A) --- -- alpha [rad]
stroke displacement
.8 --- translational velocity
***"" rotational velocity
0.,6 .. ---------- ---





S-0.2 - - Forward stroke x Backward stroke -
r -0.4 -- --
-0.6 -
-0 .8 --- -- L - I - -,- -
-1 --^ -- - T -- -- -- -'' - '-
-1.2
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
tin =2i





2 -
C3) .. .









0 0
0 -


B)o l ---g.


0.4 0.5
time/T


Figure 4-15. Kinematics of "water treading" hovering mode and the aerodynamic forces
for one complete stroke. ha=1.0c, ao=0, aa=42, k=0.5 and Re=1,700. A)
kinematics. B) lift and drag coefficients. The non-dimensional time instants
marked on the figure correspond to: A=0.03, B=0.19, C=0.31, D=0.46,
E=0.55, F=0.67, G=0.79 and H=0.95.


After reaching the maximum translational velocity and pitching angle at time/T =


0.25 (see Fig. 4-15(A)), the airfoil starts decelerating and pitching down, and a counter-


clockwise leading edge vortex starts to form (Fig. 4-16(A), at time C). The pressure drop


inside the vortex core increases (Fig. 4-16(B) at time C): therefore the lift reaches the


maximum peak during the half-stroke (- time/T = 0.3). Once the leading edge vortex is