UFDC Home  myUFDC Home  Help 



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 thermofluids 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 quasisteady models of insect flight...............................18 1.4 .2 W agner E effect .................................................................. .............. ..20 1.4.3 Leadingedge Vortex and Delayed Stall................. ............................21 1.4.4 Rotational Forces ...................................... .............. ................24 1.4.5 W ingw ake Interaction ................................................... .................. 25 1.4.6 Clap and Fling M echanism ..................................... ......... ............... 27 1.5 Objective and O utline ...................................................................... 28 2. COMPUTATIONAL METHODS FOR COUPLED FLUIDSTRUCTURE 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 asterslave C concept ........................................................... .................... 35 2.5.2 Grid Rearrangement Algorithm with Weighting Functions..................39 2.5.3 Grid Rearrangement Algorithm Based on Spring Analogy Concept........41 2.6 M em brane Structural Solver........................................... .......................... 42 2.6.1 M ooneyRivilin 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 2D 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 11. Variety of beating frequency and Reynolds number for birds and insects. ............7 12. Experimental robotic and mechanic devices for insect flapping wing study..........14 13. Principal numerical studies of insect flapping flight................... ...... .. .........18 21. Moving boundaries techniques. Advantages and disadvantages ..........................42 31. Aerodynamic performances for flexible and rigid wing. ......................................69 41. Numerical and analytical friction coefficient values for flat plate........................75 42. Grid and time increment details for the "water treading" hovering cases .............84 43. Grid and time increment details for the "normal" hovering cases ........................91 44. Kinematics parameters for the experimental and computational setup...................94 45. Kinematic parameters for "water treading" and "normal" hovering modes. The Reynolds number for both cases is 100.............. ............ ............................... 97 46. Parameters for "water treading" hovering mode employed for Reynolds number effect stu dy ..................................................... ................. 10 1 LIST OF FIGURES Figure p 11. W ing beating m otions.. .............................. ... ......................................... 12. H covering flight. .......... ... .... .. .. ..... ..... ...... .... ...... .......... .... .... 13. Wing beat frequency versus the body mass for different insect species ................9 14. Crosssectional corrugations of dragonfly wings .................................................12 15. Leading edge vortex on a hawkmoth wing. .................................. ...............15 16. Flow around a thin airfoil...............................................................................22 17. Leading edge vortex (LEV) structures at different Reynolds numbers...................24 18. Momentum transfer due to wingwake interaction..............................................26 19. Schematic of clap and fling mechanism................... ............. ............... 27 21. Illustration of the moving grid technique .............. ........ .................38 22. Masterslave algorithm performance for large rotations................................. 39 23. Current node neighbors position for grid rearrangement algorithm ....................40 24. Grid quality for the grid rearrangement algorithm............... ....................... 41 25. Generic algorithm for fluidstructure interaction .......................................... 49 31. Fixed and flapping wing M AVs. ........................................ ......................... 53 32. M A V 's w ing shapes. .............................................................................. .............54 33. Computational domain for the MAV wing. ................................. .................55 34. Streamlines, vorticity and pressure on the 6" MAV wing's upper surface as a function of the angle of attack. ...................................................... ..................57 35. W ing shape geom etry ...................................................................... .................. 59 36. Vorticity magnitude contours behind the wing at 6 angle of attack. ...................60 37. Pressure contours on the wing's upper surface and streamlines slightly above the w ing at 6 angle of attack.. .................................... .......................... .....................61 38. Pressure contours on the wing's lower surface and streamlines slightly below the w ing at 6 angle of attack. .............................................. ............................. 62 39. Pressure coefficient on the wing surface at x/c = 0.34 and 6 angle of attack. .......63 310. Pressure coefficient on the wing surface at x/c = 0.53 and 6 angle of attack........63 311. Spanwise lift distribution at 6 angle of attack. .............................................. 64 312. Spanwise drag distribution at 6 angle of attack. ................................................65 313. Surface grid for the 5" M AV wing............................................................ ........ 66 314. History of the vertical displacements of two points, one near the wing tip and the other one between the two battens. .............................................. ............... 67 315. Vertical displacement contours at various times .................................... .........68 316. Aerodynamic forces history for the 5" flexible wing MAV. ................................68 41. Numerical and analytical velocity distribution in the boundary layer along different locations on a flat plate and different Reynolds numbers. ......................................74 42. Numerical and analytical friction coefficient Cffor a flat plate versus Reynolds num ber ......................... ......................... ........................ 76 43. Skinfriction and wall velocity for an oscillating plate. ........................................77 44. Velocity distribution above an oscillating plate............... ............................... 78 45. Hummingbird in hovering flight. ........................................ ........................ 79 46. The schem atics of tw o hovering styles .............. ............. ................. ............. 80 47. Schematics of the "water treading" hovering mode ............... ..... ............... 81 48. Schematics of "normal" hovering mode. .................................. ....... ............ 82 49. The coarse multiblock grid used in the current work. The grid has 90 points around the airfoil and 70 points in the radial direction. ........................................83 410. Lift coefficient for three periods and different grid sizes and 6t=0.01....................85 411. Drag coefficient for three periods and different grid sizes and 6t=0.01..................85 412. Lift coefficient for three periods and different time increments ...........................86 413. Drag coefficient for three periods and different time increments ..........................86 414. Lift and Drag coefficients for one period ..................................... ............... 88 415. Kinematics of "water treading" hovering mode and the aerodynamic forces for one com p lete strok e................................................. ................ 89 416. Flow field snapshots of the "water treading" hovering mode...............................90 417. Lift coefficient of the "normal" hovering mode for three periods and different grid sizes and St=0.01. ......................................................................92 418. Lift coefficient history for the "normal" hovering mode. .....................................95 419. Vorticity plot for "normal hovering" mode with h/c=2.4 and Re= 115...............96 420. One cycle force history for two hovering modes.. ............................................... 98 421. Vorticity contours for two hovering modes.. ...................................................99 422. Lift coefficient for the "water treading" mode ............................................... 102 423. Vorticity contours for the "water treading" mode...................................103 424. 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 liftto 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, pressurebased NavierStokes solver along with moving grid algorithms is employed to simulate the flow field. The coupled fluidstructural dynamics of the flexible wing is treated using a hyperelastic finite element structural model, the abovementioned 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 fluidstructure interactions for flexible structures have been investigated. In the Reynolds numbers range of 7x104 to 9x104, the flexible wing exhibits selfinitiated vibrations even in steady freestream, 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 delayedstall and rapid pitchup 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 wakecapturing, 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 manmade 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 aerohydrodynamic properties and use energy more efficiently compared with those of the manmade 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 SR71 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 multidisciplinary effort involving aerohydrodynamics, advanced materials, aero/hydro/servo/elasticity, microelectromechanical 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 (SchmidtNielsen, 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. 11): 1. an outofplane motion called "flapping" (flapping angle w ); 2. an inplane 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 11. 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. 12(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 12(B). Takeoff 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 takeoff, whereas large heavy birds take long runs against the wind to reach takeoff 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 12. 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 11), 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 (11) 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, (12) 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 (13) 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 (14) Table 11. 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.5x108 370 2.0x101 (Encarsia formosa) Fruit fly 2.0x106 240 2.0x102 (Drosophila virlis) Fly 1.2x10' 190 6.4x102 (Musca domestic) Mosquito 3.5x106 320 1.3x103 (Aedes nearcticus) Honeybee 7.8x105 250 2.3x103 (Apis mellifica L.) Bumblebee Bumbutlebee 8.8x104 156 4.0x103 (Bombus terrestris) Dragonfly 7.9x104 29 4.9x103 (Anax parthenope) Locust Locust .2.0x103 20 1.0x104 (Schistocerca gregaria L.) Hummingbird 2 4 Hummingbird 2.0x102 15 1.5x104 (Patagona gigas) Sparrow 3.0x102 13 1.0x105 Pigeon 3.5x101 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 (15) 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 (16) 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 (17) 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, (18) 2 where PA is the aerodynamic power and the geometric similarity for wingspan b was used (b oc m13) (Azuma, 1992). Combining equations (17) and (18), the lower bound for wingbeat frequency is given as: f OcUm 13 o m1/6 (19) 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 (110) 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 (110) leads to: f oc m1/6, which is consistent with the result given in equation (19). 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 m1/2 and m1/6 line, as can be seen in Fig. 13. 4Wing beat frequency vs body mass for different insect species 10 , ,,,,ai different insect species I ....... .. . f = 1/2 441+ 1 Ts 1 44f ff1 4 3I 1 1 1 1 11ll 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 103 10 l Lij ^ J _~ 10 fii l6_II11__Illl_.. Figure 13. Wing beat frequency versus the body mass for different insect species. (Data collected from Azuma, 1992). I I 4IT  I1I4I1   I744 I  F F1T IT I T IT I  I4 I T II 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 FFIi 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 III 4II II4 4III I4II I I I I I I I TITI II III 108 10 7 10T6 105 10 4 10 3 102 Figure 13. 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 welldefined crosssectional 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. 14), 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 highspeed 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 14. Crosssectional 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 12. 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 largescale motion. Liu and Kawachi (1998) were the first to attempt a full unsteady NavierStokes 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 12. Experimental robotic and mechanic devices for insect flapping wing study. Author, Year Wing model Main features based on: Bennett (1970) Beetle 3D mechanical model (Melolontha flow visualization vulgaris) only wing incidence can modified do not allow force measurements Spedding and Generic 2D 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 3D mechanical model (1987), (1988) and flow visualization Savage et al. (1979) Ellington et al. (1996) Hawkmoth 3D robotic model good flow visualization (model wing =10 times insect wing) computercontrolled wing kinematic do not allow force measurements Dickinson et al. (1999) Fruit fly 3D robotic model and Lehmann (2000) good flow visualization (scaledup wing) computercontrolled wing kinematic allow measurements of shear forces using a 2D force transducer To facilitate the resolution of viscous flow at the leading edge, trailing edge and wing tip, an 00 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 spanwise flow in stabilizing the spiral leadingedge vortex as a liftenhancement mechanism, as one can see in Fig. 15. 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 remeshing 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 15. 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 secondorder 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 noninertial, 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 pitchingup 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 threedimensional, unsteady, incompressible NavierStokes 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, secondorder finite difference formulas are used on a hybrid staggered/nonstaggered grid layout. The discrete equations are integrated in time using a secondorder, dualtimestepping, 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 OH 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 twodimensional simulation of a hovering dragonfly. The solver employed a secondorder accurate central difference scheme for spatial discretization and a mixed explicitimplicit 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 3D effects, comparison of experiments and computations in 2D have also provided important insight. Wang (2000a, 2000b) demonstrated that the force dynamics of 2D 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 13. Table 13. Principal numerical studies of insect flapping flight Author, Year Fluid governing equations Method Moving boundary tracking Smith, (1996) 3D potential flow Panel method Time depended boundary conditions Liu and Kawachi, 3D NavierStokes in Artificial Grid remeshing. (1998) bodyfitted coordinates compressibility Analytical based on the initial (Chorin, 1968) grid and wing kinematics Wang, (2000a, 2000b) 2D NavierStokes for Finite difference Time depending boundary vorticity in elliptic conditions coordinates. Mittal et al., (2002) 2D NavierStokes Finite volume Immersed Boundary Method equations on fixed Cartesian grid Ramamurti and 3D NavierStokes in Finite Element Grid remeshing Sandberg, (2002) ALE formulation Spring analogy and smoothing (Arbitrary Lagrangian on unstructured grid Eulerian) Sun and Tang, (2002a, 3D NavierStokes in Artificial Time dependent coordinate 2002b) bodyfitted coordinates, compressibility transformation (Rogers et al., 1991; Rogers and Kwak, 1990) Sun and Lan, (21 114) 3D NavierStokes in Artificial Overset moving grid bodyfitted coordinates, compressibility (Rogers et al., 1991; Rogers and Kwak, 1990) Gilmanov and 3D NavierStokes Finite difference Immersed Boundary Method Sotiropoulos, (2005) on fixed Cartesian grid Mittal et al., (2006) 3D NavierStokes Finite difference Immersed Boundary Method on fixed Cartesian grid 1.4 Unsteady mechanisms identified to enhance lift in insects 1.4.1 Analytical and quasisteady models of insect flight Early models of insect flight analyzed the farfield 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 (spanwise) 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 quasisteady 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 quasisteady 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 steadystate 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 halfstroke 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 quasisteady models. Studies done on 2D 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 Leadingedge 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. 16. For a 2D 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 counterrotating 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 16. 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 2D models, the leading edge vortex in 3D 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 3D leading edge vortex stability are still unclear. Ellington and coworkers observed a steady spanwise flow from the wing's hinge to approximately threequarters 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 backsweep delta wings. However, the importance of axial flow in the stability of the leading edge vortex was questioned by experiments using a 3D 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 spanwise flow within the leading edge vortex core was quite small (25% 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 middownstroke in hovering, Liu et al. (2005) found that the leadingedge vortices show significant variations in their structures as the Reynolds number changes. Figure 17 shows the streamline structures at three Reynolds numbers, corresponding to a hawkmoth in Fig. 17(A), a fruit fly in Fig. 17(B), and a thrips in Fig. 17(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 leadingedge vortex is observed on the paired wings with a steady spanwise flow at the vortex core, breaking down at approximately threequarters of the span towards the tip. At a Reynolds number of 120, corresponding to a fly (Fig. 17(B)), the breakdown apparently disappears, and a leadingedge vortex is found to be connected to the tip vortex, forming a longer, leading edge tip vortex. The vortex shows strong stability, while the spanwise 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 leadingedge vortex, the tip vortex and the trailing vortex is observed (Fig. 17(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 pressuregradient, the centrifugal force, and the coriolis forces are all together responsible for the leadingedge vortex stability. Their roles at different Reynolds numbers should be studied to help shed light on these findings. Figure 17. 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 spanwise 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 spanwise axis near the end of every stroke in order to obtain a positive angle of attack for the next halfstroke 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 halfstroke 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 Wingwake 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. 18. 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 pitchingup rotation of the wing, and not to the wingwake 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 18. Momentum transfer due to wingwake 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 takeoff. However, the contribution of the wingwake 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 clapandfling mechanism is the physical interaction of the left and right wing during dorsal stroke reversal. This mechanism was proposed by WeisFough (1973) to explain the lift generation in the chalcid wasp. Its schematic is presented in Fig. 19. clap a b c d fling v a b c d Figure 19. Schematic of clap and fling mechanism. Based on WeisFough (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 buildup of circulation. A detailed theoretical analysis of the clapandfling 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 (10104), with focus on twodimensional 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 manmade 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 fluidstructure 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 pressurebased fluid solver. The finite element structural solver and the governing equations for solids are also presented. The fluidstructure 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 fluidstructure interaction problem is presented. In chapter 3, the flow over fixed rigid and flexible wing MAVs is analyzed. The low aspectratio 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 twodimensional 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 delayedstall and wakecapturing mechanisms. Finally we conclude with the summary and the directions for the future work in chapter 5. CHAPTER 2 COMPUTATIONAL METHODS FOR COUPLED FLUIDSTRUCTURE INTERACTIONS 2.1 Overview This study investigates the aerodynamics of fixedwing 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 wingwake 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 PressureImplicit 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 threedimensional NavierStokes equations for incompressible flow written in curvilinear coordinates are solved. Two popular methods to solve the NavierStokes 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 pressurebased 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 remesh 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 threedimensional NavierStokes equations for incompressible flow in Cartesian coordinates can be written, using indicial notation, as follows: ai+ a ( )=O (21) at cx Ca(u)+ I + u,) (22) 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 (23) Oxi cx 3 Ox, where / is the molecular viscosity. For arbitrary shaped geometries, the NavierStokes 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 (24) 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 (25) _z z z _ To solve for the NavierStokes 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 pressurecorrection equation. The expressions for the metrics, the determinant of the Jacobian transformation matrix, and the fluxes at the cell faces for the three dimensional NavierStokes equations in the generalized bodyfitted 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 bodyfitted 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 nonstaggered 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 divergencefree velocity field within a desired convergence tolerance, therefore enforcing the pressurevelocity 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 noniterative process. Here the pressurevelocity 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 NavierStokes equations, proper boundary conditions are required. At the interface between the fluid and solid structure, the noslip 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 bodyfitted 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 (25). 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 (26) 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 (27) Ot 9^ t0 B7 Wt 8 OR Integrating equation (27) with a firstorder, 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 Masterslave 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 remeshing 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 remesh the multiblock structured grid for fluidstructure interaction problems. For multiblock structured grids, for simplicity, CFD solvers often require pointmatched grid block interfaces. This method is based on the masterslave concept and maintains a pointmatched grid block interface while maintaining grid quality and preventing potential grid crossover. 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 offbody points. One difficulty for a multiblock grid resides in the way in which the vertices of each block are moved, if a pointtopoint 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. 21(A), the edge and interior points can be obtained by a 3stage 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. 21(A), the interpolation can cause a discontinuity at the interface. To avoid creating undesirable grid discontinuities, the offbody vertices of a grid block are linked to a surface point and thus move in a similar way. Therefore, for each offbody 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 (28) 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) (29) 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, ) RUMI"'C" " = i:i: I% ' 7:  .7 :i.i..iiii:.i Figure 21. 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} (210) 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. (211) 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. 21(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 22. Masterslave algorithm performance for large rotations The masterslave 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. 22, for large deformations, the grid quality deteriorates, leading to negative cell volumes. 2.5.2 Grid Rearrangement 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. 23), / 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 23. Current node neighbors position for grid rearrangement algorithm. 41 vi) at block interfaces the displacement is interpolated as the average between neighboring blocks values, in order to maintain the pointtopoint 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. 24(A). 14 i 12 "10 " 12 1 14 0 1 15 10 5 0 10 I 1 A) x B) x Figure 24. Grid quality for the grid rearrangement algorithm. A) Close to the body. B) Overall domain. The disadvantages of this algorithm are the large skewness near the outer boundary (Fig. 24(B)) imposed by the nonmovement 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 Rearrangement Algorithm Based on Spring Analogy Concept To overcome the main disadvantage of the grid rearrangement algorithm described in section 2.5.2 and of the masterslave algorithm (section 2.5.1), a new algorithm based on the springanalogy 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 crossover encountered in the original masterslave algorithm (Fig. 22) is eliminated. Table 21. Moving boundaries techniques. Advantages and disadvantages. No. Grid Advantages Disadvantages Observations movement algorithm 1. Masterslave 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 2D 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 2D 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 21. 2.6 Membrane Structural Solver The main characteristic of membrane structures is their large deformations implying a nonlinear stressstrain 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 straindisplacement 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 Fopplvon 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 twodimensional 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, twodimensional, 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 twodimensional case, the threedimensional membrane model introduces several complicating factors. In the twodimensional "sail equation", the tension is a single constant. In threedimensional 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 rubberlike membrane. Ding et al. (2003) numerically studied partially wrinkled membranes. The threedimensional 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 MooneyRivilin model (Mooney, 1940). A brief review of their membrane model is given next. 2.6.1 MooneyRivlin Model The MooneyRivlin 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) (212) 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) (213) where cl and c2 are two material constants. A material that obeys equation (213) is known as a MooneyRivlin material. For an initially isotropic, incompressible membrane, the general stressstrain relation is written as: S= pC +2[(c, +c21)IcC] (214) where S is the second PiolaKirchhoff 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. (215) The hydrostatic pressure is determined by the condition that S33 = 0, and the formula is: p=2(c +C2I1 C2 )1 (216) where X3 is defined by = 33 (t)= h(t)/ho (217) 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 GreenLagrange 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), (218) 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 (218) 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 (219) 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 (219)), 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 WilsonO method (Wilson, 1973) is used. It is a onestep method which is easier to code and behaves well in fluidstructure interaction problems. The Wilson0 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 secondorder 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 fluidstructure interaction problems, the CFD grid needs to be automatically regenerated 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 onedimensional case can be written as: N H(x) ax x x 12logx x, (220) 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, (221) 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, Remeshing Based on the interpolated displacements update the CFD grid Recompute the CFD grid metric terms Update the Jacobian to satisfy the GCL Next time step, t = t + At STOP Figure 25. Generic algorithm for fluidstructure interaction In order to conduct timeaccurate computations for fluidmembrane dynamics, iterations need to be done at each time instant to enforce the fluid and structure coupling. A generic algorithm for fluidstructure 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 fluidstructure coupling, is shown in Fig. 25. 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 1020 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 biochemical 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 (PomsinSirirak 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 (104105) 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 microelectronicmechanical systems (MEMS) for MAVs are discussed by GadelHak (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 31(A) demonstrates a fixed wing design by Ifju and coworkers 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 chordwise 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 PornsinSirirak et al. (2000) and is shown in Figure 31(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 31. 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 32 shows the wing shapes designed by Ifju and coworkers for their MAVs and are used as a base for the present numerical simulations. Figure 32(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 32(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 crosssection and a dihedral angle of almost 7. The double curvature airfoil crosssection offers better transversal (pitch) stability characteristics for flying wing airfoils. A) I B)Z Figure 32. 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, NavierStokes equations for the incompressible flow are solved as described in Chapter 2. A pressurebased algorithm is used to solve the threedimensional NavierStokes equations written in curvilinear coordinates. Both firstorder and second order upwind schemes are employed for convection terms and secondorder 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, multiblock 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 33. (6c, 5c, 6c) J 4 seC) (6 56 OUTLET .",, ~~, ^ OUTLET . c 5c, 6c) (c A INLET (9c, 5c, 6c) B INLET (1., .,4e) Figure 33. 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 chordwise 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 zerogradient 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 highpressure region below the wing surface and the lowpressure 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 Te AR (31) 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 32, 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 lowaspect wings, it is important to investigate the tip vortex effects on the wing aerodynamics. In general, tip vortex effects are twofolded: * Tip vortex causes downwash that decreases the effective angle of attack and increases the drag force (Anderson, 1989). * Tip vortex forms a lowpressure 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 outerboundaries of the computational domain is presented in Figure 33(A). Based on the root chord length and a freestream 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 34 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 lowpressure area on the wing's upper surface becomes larger. A) B) Figure 34. 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 lowpressure 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 lowaspect 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 lifttodrag ratio. To accomplish this objective we try to maintain/increase the lift while decreasing the induceddrag. 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. 35(A)) with dimensions specified in section 3.2.1; (2) the modified wing (Fig. 35(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. 35(C). Original wing shape Modified MAV wing Modified MAV wing and Endplates A B C Figure 35. 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. 33(A). The boundary conditions are described in section 3.2.2, with the addition of noslip boundary conditions assigned to the endplate surfaces. The freestream 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 underpredicts the lifttodrag 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. 36, the vorticity magnitude contours behind the wing are plotted in planes perpendicular to the xaxis, i.e., the streamwise 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 36. 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 highpressure region on the lower wing surface from reaching the lowpressure 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 lowpressure zone present on the upper wing surface, as indicated by the streamlines plotted in Fig. 37(A), (B). The lowpressure zone associated with the tip vortex core can be clearly identified in the plot of the spanwise pressure variation on the wing's upper surface (Fig. 39(A) and Fig. 310(A)). Wing tip   : p 9079685746352413 2 9 20 B Wing tip as_ r Wng tip ._  K p: 90 7968 5746 3524 13 2 9 20 Figure 37. 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 highpressure zone (below the wing surface) and the lowpressure zone (above the wing surface) can also be observed by looking at the streamlines slightly below the wing surface, as shown in Fig. 38(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. 38(B). The decrease of the highpressure area due to the circulatory motion in the absence of the p 90 79 685746 35 24 13 2 9 20 I I _ 62 endplate can also be recognized by looking at the spanwise pressure coefficient on the lower wing surface, which is plotted in Fig. 39(B) and Fig. 310(B). __ VWng tip A p: 302418126 0 6 12 18 24 30 p 30241812 6 0 6 12 18 24 30 LW Vng tip p 302418126 0 6 12 18 24 30 SVWng tip p: 30 24 18 12 6 0 6 12 18 24 30 Figure 38. 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. 39(A) and Fig. 39(B)), reducing the vortex lift. On the other hand, a lower velocity slightly below the wing increases the highpressure 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 highpressure zone on the lower wing surface in the presence of the endplate can be clearly seen from the spanwise pressure coefficient on the lower wing 63 surface plot (Fig. 39(B) and Fig. 310(B)), and also from the pressure contours on the lower wing surface shown in Fig. 38. 02 0A) A) 01 02 03 Z / Local span 0a C 04 05 Z / Local span Figure 39. 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 90 05  01 015  02  01 02 03 Z / Local span 1 0 : 04 05 01 02 03 Z/ Local span Figure 310. 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 spanwise locations. The result is the lift and drag estimation in the spanwise distribution. Figure 311 shows the spanwise 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 311. Spanwise 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 spanwise drag distribution plotted in Fig. 312, 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 liftto 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 lifttodrag 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 312. Spanwise 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. 32(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 carbonfiber battens, as shown in Fig. 313(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. 313. Grey area = Rgid part, Red = Wing tip regions; Green = Rigid leading edge boundary; Black = Rigid central panel Blue = Batten locations Figure 313. 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 freestream 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 2x104 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 fluidstructure framework is presented in chapter 2. Even though the flow stream is uniform, the membrane exhibits selfinitiated vibrations, as one can see in Fig. 314. 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 o005v  near wingtip  between the battens oil 1 L 0 001 002 003 004 005 006 007 008 009 01 time [sec] Figure 314. 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. 315 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. 314. tme =0005 *,. time=015 00 le025 a J1 U:::^ Ju^ 4 ...... &L.. : & L i .B0045 ., 0056 Y u ? Figure 315. 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. 316 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 316. 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 31. Table 31. 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 lifttodrag 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 steadystate, laminar, incompressible NavierStokes 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 NavierStokes solver was coupled with a finite element structural solver to investigate the fluidstructure 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 freestream. For small to moderate angles of attack, and before stall limit, the timeaveraged 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 threedimensional effects are critical for MAVs, experiments and computations in twodimensions continue to provide important insight into flapping wing aerodynamics. From the computational point of view, 2D simulations require less computational power, and grid sensitivity and convergence studies can be performed more readily (compared to 3D cases). Wang (2000a, 2000b) demonstrated that the aerodynamics of 2D wings could help explain some of the enhanced lift observed in insect flight. In this chapter, the flow structures around a hovering 2D 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 (41) On) all In the present approach, equation 41 is discretized using a first order scheme: Twal = (Utan )p (tan)wall (42) 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 (nonzero 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 NavierStokes equations are solved on a grid with 100 points along the streamwise and 60 points along the crossstream direction. The distance from the wall to the first cell center is 5x104c, 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 freestream 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 noslip boundary condition is imposed on the airfoil surface. In Fig. 41, 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) (43) where Um is the freestream 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 41. 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 41 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 semiinfinite 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 trailingedge of a finitelength plate. The discrepancy between the analytical and numerical velocity distribution is visible at locations near the trailingedge of the plate, as one can observe in Fig. 41 for all Reynolds numbers. Figure 42 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 skinfriction coefficient for low Reynolds number. The analytical expressions for skinfriction coefficient, defined for 2D cases as Cf = Drag/(0.5pcU2) are presented in Table 41. Table 41. 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 42 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 I4 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 42. 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 subsection, 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 nonzero). If we consider an infinite plate oscillating along the xaxis, the wall velocity is defined as: u(y = 0, t) = Uo cos(ot) (44) where Uo is the maximum velocity and c) = 27rf ; being the oscillation frequency. The integration of equation 44 gives the plate displacement: h(t) = h, sin()t); I i it h = Uo l/ (45) The analytical velocity field for this motion is given by Stokes (1851): u U/U = exp(r) cos(ot r/) (46) where r7 = yJ) /2v is a nondimensional crossstream coordinate. Following the definition given in equation 41, the analytical wall shear stress for an oscillating plate is given by (White, 1991): Trwa = U0o p sin(t r/4) (47) 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 NavierStokes 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 43. Skinfriction 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. 43, the numerically predicted wall shear stress on the finite plate is compared with analytical values for an infinite plate given by equation 47. 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. 43 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 44. 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 46, is plotted, along with the computed velocity field, in Fig. 44 for different time instants during one halfperiod. The skin friction and velocity profiles presented in Fig. 43 and Fig. 44 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 2D 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 3D) 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 45 shows the figure "8" pattern observed in hovering hummingbirds. Figure 45. Hummingbird in hovering flight. A) Forward stroke, B) Backward stroke. Picture taken by Wei Shyy. In an extensive survey, WeisFogh (1973) noticed that most insects, such as fruit flies, bees, and beetles, employ symmetric backandforth strokes near a horizontal stroke plane, as seen in Fig. 46(A). WeisFogh 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 46. 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 socalled "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), (48) and a rotational component (pitching), given by: a(t) = ao + a, sin(2kt). (49) In equation 48, ha is the plunging amplitude and k is the reduced frequency of the sinusoidal oscillation. In equation 49, 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 (410) 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. 47. Y Forward stroke ha ha Backward stroke ha ha Figure 47. 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) (411) a(t) = ao a sin(2kt) (rotation) (412) 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 410. A schematic of the motion is shown in Fig. 48. Forward stroke X Backward stroke X ha7 ha Figure 48. 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 bodycentered 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.4x103c. In this work, an Otype 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. 49. The minimum grid spacing adjacent to the airfoil varies from 2.4x103c (coarsest grid) to 1.0x103 (finer grid). DgrO (90x70, dy=2.4e03) > o 4 10  10 0 10 x Figure 49. The coarse multiblock 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, NavierStokes 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 48 and 49 with the same kinematics parameters as in the "water treading" hovering mode presented earlier. Table 42 provides a summary of the cases run, along with grid and time increment details. A schematic of the movement is presented in Fig. 47. Table 42. 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.0x10c) 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 410 and 411 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 412 and 413 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 410. 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 411. 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 412. 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 413. 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. 410, 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 pitchdown 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 414(A) shows the lift and drag coefficients history for one period as obtained by Liu and Kawachi, while Fig. 414(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 NavierStokes 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 1x102 is used for this analysis. Eight time instants are selected for analysis during one hovering period and are identified in Fig. 415. The airfoil position and pitching angle, as well as the translational and rotational velocities, are plotted in Fig. 415(A). Fig. 415(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 414. 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 pitchup rapidly. Fig.415(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. 415(A)). Figure 415(B) shows that, during this period, the lift curve reaches a plateau. The clockwise vortex shed earlier (Fig. 416(A) at time A) accelerates the flow over the upper surface of the airfoil, and creates a lowpressure zone (Fig. 416(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 ..   S0.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 415. 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 nondimensional 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. 415(A)), the airfoil starts decelerating and pitching down, and a counter clockwise leading edge vortex starts to form (Fig. 416(A), at time C). The pressure drop inside the vortex core increases (Fig. 416(B) at time C): therefore the lift reaches the maximum peak during the halfstroke ( time/T = 0.3). Once the leading edge vortex is 