When vortices stick: an aerodynamic transition in tiny insect flight

Insect flight has been a topic of considerable interest in a variety of fields. For instance, evolutionary biologists have long been fascinated by the evolution of flight and the functional morphology of flight apparatus(Dudley, 2000). Biologists have also been interested in the diversity of flight kinematics, flight control and aerodynamic mechanisms. More recently, engineers have been involved in the study of insect flight and its aerodynamic mechanisms as interest has emerged in the design of micro-robotic flyers(Ellington, 1999). There has also been progress in mathematics and computer science to understand better the complicated fluid dynamics involved in insect flight, with use of analytical models and computational fluid dynamics (CFD; Lighthill, 1973; Sun and Tang, 2002; Wang, 2000a,b). In early studies, scientists attempted to describe the lift-generating mechanisms of insect flight with traditional quasi-steady-state aerodynamic theory (Ellington, 1984a; Lighthill, 1975; Osborne, 1951). This method essentially reduces the motion of the wing to a series of consecutive states of steady flow over the entire wing or span-wise sections of the wing. However, early quasi-steady-state analysis often failed to predict the magnitude and direction of forces measured directly in tethered insects(Cloupeau et al., 1979; Ellington, 1995; Zanker and Götz, 1990). These findings have led researchers to explore further possible unsteady mechanisms of lift and thrust generation and to develop revised quasi-steady-state models (Sane and Dickinson, 2002).

One unsteady mechanism that has been suggested as a means to provide additional lift during insect flight is delayed stall. Essentially, a large leading edge vortex (LEV) is formed at the beginning of each half stroke and remains attached to the wing until the beginning of the next half stroke. In a typical airplane wing translating at a constant angle of attack, stall occurs above some critical angle of attack when the LEV is shed and lift forces consequently drop. Stall, however, appears to be delayed or suppressed for revolving insect wings operating at high angles of attack. The question then remains as to whether or not the LEV would be shed during translation at some time beyond the length of the half stroke or whether it would remain attached indefinitely in a three-dimensional (3-D) stroke. Attached LEVs have been observed in flow visualization studies of the hawkmoth, a dynamically scaled model flapper (Ellington et al.,1996; van den Berg and Ellington, 1997a,b; Willmott et al., 1997), and revolving model wings (Usherwood and Ellington, 2002). An attached LEV was observed in a 3-D CFD simulation of hawkmoth flight by Liu et al.(1998), and Birch and Dickinson (2003) also observed a stable attached LEV using time-resolved digital particle image velocimetry(DPIV) of the flow around the wings of a dynamically scaled robotic insect. In a later study, Birch et al.(2004) showed that this stable attached LEV is a robust phenomenon for Reynolds numbers (Re) in the range of 120 to 1400.

To understand the mechanism of the formation and attachment of the LEV,consider a wing translated from rest immersed in a viscous fluid initially at rest. At the onset of motion, the solid body has a non-zero tangential velocity relative to the surrounding fluid. This motion shears the fluid, and the discontinuity creates a sheet of concentrated vorticity. At later times,the vortex sheet is transported away from the boundary by diffusion and convection. As a result of the negative pressure region generated instantaneously behind the wing due to the motion of the fluid, the vortex sheet `rolls up' to form the LEV. The `attachment' of the LEV to the wing maintains the negative pressure region behind the wing, which leads to higher lift forces. A question of interest to both fluid dynamists and biologists is why the shedding of the LEV is delayed or suppressed at high angles of attack during insect flight. Several authors have suggested that axial flow along the wing derived from a span-wise pressure gradient stabilizes the LEV and delays stall (Ellington, 1999; Liu et al., 1998).

Although there is no rigorous theory regarding the stability of the attached LEV, insight can be gained by considering a 3-D fixed wing in a steady flow with an attached LEV. There are three processes occurring in this case: the convection of the vortex, the intensification of vorticity when vortex lines are stretched, and the diffusion of vorticity by viscosity(Acheson, 1990). In order for the LEV to be steady and remain attached to the wing, these three processes should be balanced. Note that in the two-dimensional (2-D) case, vortex stretching cannot occur. This might account for differences observed in the stability of the leading edge vortex between the 2-D and 3-D cases. In addition to differences in the dimensionality of the problem, Birch et al.(2004) found that the structure of the stable attached LEV differed for high (1400) and low (120) Re. At an Re of 1400, they found an intense narrow region of span-wise flow within the LEV. At an Re of 120, this region of span-wise flow was absent, suggesting that the 3-D mechanism contributing to the stability of the LEV takes different forms at high and low Re.

Weis-Fogh (1973) proposed another unsteady lift-generating mechanism termed `clap-and-fling', which is mostly observed in the smallest flying insects(Ellington, 1984b; Weis-Fogh, 1975). Basically,the wings clap together at the end of the upstroke and are then quickly peeled apart at the beginning of the downstroke. This motion has been shown both experimentally and analytically to enhance circulation around the wings and augment the lift generated during the downstroke(Lighthill, 1973; Spedding and Maxworthy,1986).

Another possible mechanism for enhanced lift generation in insect flight is that circulation around the wing is enhanced by the quick rotation of the wing at the end of the downstroke. Dickinson et al.(1999) suggested that large rotational forces generated during rotation induce a net lift force that is analogous to the Magnus effect seen in the case of a spinning baseball. In this case, however, the force is generated by the rotation of a flat plate rather than a round cylinder, and the net rotational force acts approximately normal to the chord of the wing. Sane and Dickinson(2002) incorporated this idea(based on thin airfoil theory) into a quasi-steady-state model of flapping flight using rotational force measurements taken from a dynamically scaled model insect. Sun and Tang(2002) suggest that these large rotational forces can be attributed to the rapid generation of vorticity during wing rotation and reversal. Walker(2002) suggests that the large forces generated during wing rotation can be described by a quasi-steady-state model without a rotational term analogous to the Magnus effect.

`Wake capture' and vortex effects from previous strokes are other possible mechanisms proposed to generate lift during insect flight(Dickinson, 1994; Sane and Dickinson, 2002; Wang, 2000b). Essentially,vortices produced from previous strokes enhance the lift generated by subsequent strokes. One way this might act to enhance lift is that the flow generated by one stroke increases the effective fluid velocity at the start of the next stroke. By definition, these forces are not observed during the first stroke. As one would expect, they depend upon the point of rotation, timing of rotation and rotational speed (Dickinson,1994; Dickinson et al.,1999; Wang,2000b). Wang(2000a,b)described this phenomenon computationally and found that there exists an optimal flapping frequency for lift generation. This optimum results from two time scales: vortex growth and the shedding of the LEV. In a 3-D numerical simulation, Sun and Tang(2002) did not find evidence for lift augmentation via wake capture and argue that enhanced lift attributed to wake capture can be explained by inertial forces. However, Birch and Dickinson (2003) showed experimentally that wake capture can influence lift forces based on the magnitude and distribution of vorticity during stroke reversal.

Although there have been a number of theoretical and experimental studies investigating lift generation in larger insects, few have considered those insects that fly at Re below 100. These insects are therefore said to be in the `twilight zone' of flight(Dudley, 2000). The lack of emphasis on these small insects could be partially due to several experimental and mathematical difficulties. For example, it would be rather difficult to measure actual lift and drag forces for insects this small. Kinematic analyses using video are expensive given the high range of wingbeat frequencies estimated for tiny insects. For example, measured wingbeat frequencies can be as high as 1046 Hz (Sotavalta,1947). Furthermore, most analytical work assumes that the fluid is inviscid and it seems unlikely that this is a good approximation for Re in the range of 10 to 100. Experimentalists have proposed several ideas as to how these insects generate lift. One idea involves an asymmetric stroke using a mechanism similar to that which generates thrust in rowing. Lift is generated during the downstroke, and the wing is turned to minimize negative lift on the upstroke (Horridge,1956; Thompson,1917). Another idea is that lift enhancement from clap-and-fling is sufficient for flight in this regime(Weis-Fogh, 1973).

There are several reasons to believe that flight aerodynamics change significantly for Re below 100. It is well known that below an Re of ∼40, vortices are no longer shed behind cylinders immersed in a moving fluid (Landau and Lifshitz,1959). Experimental work shows that this is also true for fixed plates (Batchelor, 1967). This transition might alter the lift enhancement generated by wake capture. In fact, Wang (2000b) found that the lift coefficients are more than halved for flapping 2-D wings when the Re is lowered from 157 to 15.7, and Wu and Sun(2004) found in a 3-D simulation that lift coefficients were decreased while drag coefficients were greatly increased for Re below 100. Walker(2002) also argues that for low Re flight, viscous forces become increasingly important to the force balance. Furthermore, at some critical Re, separation at the leading edge of the wing does not occur and the LEV does not form. The authors assert that, in addition, the trailing edge vortex will not form below this critical Re.

In the present study, the immersed boundary method was used to simulate a simple two-dimensional wing during one complete stroke for Re ranging from 8 to 128. These simulations were constructed to be similar to Dickinson and Götz's experiments (Dickinson,1994; Dickinson and Götz,1993) using a single dynamically scaled robotic wing. The motion of this wing was strictly 2-D and was divided into three separate stages:translation, rotation and translation back through the previous stroke. However, later experiments by Dickinson et al.(1999) on a fully 3-D robofly did not separate the rotation of the wing from the translational phase. The motion used in our simulations is a 2-D version of that used in the later experiments. Consequently, it is also similar to the motion used by Sun and Tang (2002) and Ramamurti and Sandberg (2002), who modeled Dickinson's experiments with CFD. The purpose of our simulations, however, is not to mimic previous work but rather to investigate what happens when the Re is lowered to that of the smallest flying insects using the same wing kinematics. We do, however, compare results with published lift and drag data for Re ranging from ∼6 to 200.

where ρ is the density of the fluid, l is a characteristic length of the wing, U is the velocity of the fluid, μ is the dynamic viscosity and ν is the kinematic viscosity of the fluid. Dickinson and Götz used an aluminum wing with a chord of 5 cm immersed in a sucrose solution with a kinematic viscosity of 0.0000235 kg m^-1s^-1 (∼20 times that of water) moving with a characteristic velocity in the range of 0.04-0.12 m s^-1. In addition, the dimensions of the sucrose tank used in the physical experiment were 1 m in length by 0.4 m in width. The same parameters were used in all of the following numerical experiments with two exceptions: (1) the size of the computational grid was increased to 1 m×1 m to reduce wall effects at lower Re and (2) the translational velocity was changed to vary the Re.

The motion of the model wing is a simplification of flight in D. melanogaster. The `downstroke' is defined as the motion of the wing from the dorsal to the ventral side of the body, and the `upstroke' is the motion from ventral to dorsal (see Fig. 1). The body of the insect is tilted upright during flight so that the flapping motion of the wing is approximately horizontal. The motion of the downstroke is divided into three stages: (1) translational acceleration at the beginning of the downstroke, (2) constant translational velocity and constant angle of attack during the middle of the downstroke and (3) translational deceleration and rotation at the end of the downstroke. Similarly, the motion of the upstroke is divided into three stages: (1) translational acceleration and the end of rotation at the beginning of the upstroke, (2) constant translational velocity and constant angle of attack during the middle of the upstroke and (3) translational deceleration at the end of the upstroke. Throughout the paper, `stroke' is defined as an entire stroke cycle. `Half stroke' refers to one downstroke or one upstroke (half of the entire stroke cycle). In all simulations, the center of rotation is located 0.2 chord lengths from the leading edge of the wing.

where τ_d is the dimensionless time when deceleration begins,and τ_final is the dimensionless duration of the entire stroke. The translational velocity during the upstroke is symmetric to the downstroke and may be constructed similarly. Unless otherwise noted, τ_dwas taken to be 10.8 (this gives a translation of ∼5 chords during both the downstroke and upstroke, which is similar to that occurring in Drosophila flight), Δτ_trans was taken to be 0.65, and V ranged from ∼0.00375 to 0.06 m s^-1.

The immersed boundary method (Peskin,2002) was used to calculate the flow around the wing. The essence of this method is that the deformation of a flexible boundary generates forces on the fluid, and the boundary itself moves at the local fluid velocity. For these simulations, we wanted the wing to move with small deformations in a prescribed motion. To achieve this, a target boundary that does not interact with the fluid is attached with virtual springs to the actual boundary. This target boundary moves with the desired motion, and the spring attachments apply a force to the actual boundary that is proportional to the distance between the two (see Fig. 3). The force applied to the boundary by the target boundary and the deformation of the boundary are then used to calculate the force applied to the fluid.

where u(x, t) is the fluid velocity, p(x, t)is the pressure, F(x, t) is the force per unit volume applied to the fluid by the immersed wing, ρ is the density of the fluid,and μ is the dynamic viscosity of the fluid. The independent variables are the time t and the position x. Note that bold letters represent vector quantities. Equations 8 and 9 are the Navier-Stokes equations for viscous flow in Eulerian form. Equation 9 is the condition that the fluid is incompressible.

where f(r, t) is the force per unit area applied by the wing to the fluid as a function of Lagrangian position and time,δ(x) is a two-dimensional delta function, and X(r,t) gives the Cartesian coordinates at time t of the material point labeled by the Lagrangian parameter r. Equation 10 spreads force from the boundary to the fluid grid, and equation 11 interpolates the local fluid velocity at the boundary. The boundary is then moved at the local fluid velocity, and this enforces the no-slip condition. Each of these equations involves a two-dimensional Dirac delta function δ, which acts in each case as the kernel of an integral transformation. These equations convert Lagrangian variables to Eulerian variables and vice versa.

These equations describe the forces applied to the fluid by the boundary in Lagrangian coordinates. Equation 12 describes the force applied to the fluid as a result of the target boundary. f_targ(r, t) is the force per unit area, k_targ is a stiffness coefficient, c_targ is a damping coefficient, and Y(r,t) is the position of the target boundary. Equation 13 describes the force applied to the fluid as a result of the deformation of the actual boundary, which is here modeled as a beam. f_beam(r,t) is the force per unit area, and k_beam is a stiffness coefficient. Equation 14describes the force applied to the fluid as a result of the resistance to stretching by boundary given as f_str(r, t), where k_str is the corresponding stiffness coefficient. Finally, Equation 15 describes the total force applied to the fluid per unit area, f(r, t), as a result of both the target boundary and the deformation of the boundary.

The discretization of the immersed boundary method used in these simulations has been described before in depth(Peskin and McQueen, 1996). We did, however, make one change to the method described in that paper. The operator

where C_L is the lift coefficient, C_D is the drag coefficient, S is the surface area per unit length of the model wing, U is the velocity of the boundary, F_D is the drag force per unit length, F_L is the lift force per unit length, and ρ is the density of the fluid. In the 2-D case, the surface area of the boundary means the area of a rectangle with width equal to the chord length of the wing and length equal to the unit distance (in this case, 1 m). Therefore, Sis just the chord length of the wing. It should be noted that these definitions were derived for high Re flows. For Re well below 1, force scales as μlU, where l is some characteristic length, μ is dynamic viscosity and U is velocity. For intermediate Re, forces on the boundary scale as some combination of the high and low Re approximations. However, we use the high Re convention for comparison with other results and note that C_D and C_L become functions of Re as the Re decreases.

We considered one set of stroke kinematics and varied the Re by changing the speed of translation and rotation of the wing. For these simulations, the angle θ was set to 135° during the downstroke and 45° during the upstroke, so that the angle of attack, α, would be 45° in both cases. Since the stroke is symmetric, the downstroke may also be thought of as the first half stroke, and the upstroke may be thought of as the second. Each simulation considered only the first stroke cycle. Therefore,the steady state of 2-D flapping flight will differ from the results of these simulations.

Drag coefficients as functions of time (expressed as the fraction of the stroke) for Re ranging from 8 to 128 are plotted in Fig. 4. The drag coefficient increases as Re decreases. This dependence on Re is expected in this intermediate range since the high Re approximation for the drag force no longer applies. The drag coefficient peaks during the advanced rotation of the wing. It also increases during times when the wing is accelerated. These variations in drag coefficient are consistent with the experimental results of Dickinson et al.(1999) and the computational results of Sun and Tang(2002). Fig. 5 shows the drag coefficients averaged during pure translation for the downstroke and upstroke for Re ranging from 8 to 128. For comparison with experimental results, steady-state drag coefficients measured by Thom and Swart(1940), mean drag coefficients of a wing translated from rest measured by Dickinson and Götz(1993), and mean drag coefficients of a wing translating through its wake measured by Dickinson(1994) are also plotted. Mean drag coefficients are larger during the upstroke for each Re. This phenomenon could be explained by the fact that the wing travels through its wake during the upstroke, increasing the velocity of the wing relative to the fluid. Similar results were found by Birch and Dickinson(2003) when drag coefficients were compared for the first and second half strokes using a dynamically scaled robotic insect. They found that the velocity of the fluid relative to the wing was greater at the beginning of the half stroke as the wing travels through its wake, resulting in larger drag forces.

Lift coefficients as functions of time (fraction of the stroke) are plotted in Fig. 6. The variations in lift for the different Re can be divided into two groups. For Re of ≥64, lift peaks during the initial acceleration of the wing. During pure translation for the downstroke, lift coefficients begin to oscillate. Large lift coefficients are generated during wing rotation and the subsequent acceleration of the wing at the beginning of the upstroke. During the upstroke translation, stronger oscillations in the lift coefficients are shown. For Re of ≤32, lift coefficients peak during acceleration and drop to relatively constant values during pure translation. Lift coefficients peak again during the rotation of the wing and subsequent acceleration at the beginning of the upstroke. After acceleration, the lift coefficients then drop to relatively constant values during the pure translation phase of the upstroke. Small oscillations, however, begin to grow in the Re∼32 case. This Re appears to be on the border of a transition that will be discussed in more detail below. What may not be apparent from this plot is that lift coefficients would continue to oscillate for Re of ≥64 if translation continued, but lift coefficients for Re of ≤32 would settle to constant values. Increased lift during acceleration and rotation is consistent with the results of Dickinson et al.(1999) and Sun and Tang(2002).

Mean and peak lift coefficients during downstroke translation are plotted in Fig. 7. Experimentally determined mean lift coefficients are also plotted(Dickinson and Götz, 1993; Thom and Swart, 1940). Mean lift values for Re of 8 and 16 are slightly larger than those measured by Thom and Swart. Their experimental values, however, were measured during steady translation. Downstroke lift coefficients shown in Fig. 7 seem to approach these experimental values. Fig. 8shows the average lift to drag ratio during translation for the downstroke as a function of Re. Lift/drag increases with increasing Re.

where u(x, t) and v(x, t) are components of the fluid velocity: u(x, t)=[u(x, t), v(x, t)]. The density of the streamlines is proportional to the speed of the flow.

For an Re of 128, vortex shedding plays an important role in the variation of lift throughout the stroke. In Fig. 9A,B, it is easy to see that an attached LEV has formed while the trailing edge vortex is being shed. This corresponds to a growth in lift forces. In Fig. 9C, the leading edge vortex is being shed while a new trailing edge vortex is formed. This corresponds to a drop in lift. During rotation, the leading and trailing edge vortices are shed (Fig. 9D,E). After rotation, the wing moves back through its wake(Fig. 9F-I). At the beginning of the upstroke, a new LEV is formed and a new trailing edge vortex is formed and shed (Fig. 9E,F). This corresponds to an increase in lift. In Fig. 9G, the LEV is shed and a new trailing edge vortex is formed. This results in a drop in lift. Finally, a second LEV is formed and the trailing edge vortex is shed, resulting in another lift peak(Fig. 9I). It has been shown by several studies that in actual insect flight the LEV remains attached to the wing for the duration of each half stroke, and the trailing edge vortex is shed. This sustained vortical asymmetry (attached LEV and shed trailing edge vortex) results in higher lift forces(Birch and Dickinson, 2003; Ellington et al., 1996). The fact that the LEV is not stable in our 2-D simulations supports the idea that the LEV in three dimensions is stabilized by span-wise flow.

For an Re of 8 (Fig. 10), leading and trailing edge vortices remain attached throughout the downstroke. Fig. 10A-Cshows the streamlines around the wing during the downstroke. Vortices form on the leading and trailing edges of the wing and remain attached until the end of the downstroke. Since no vortices are shed, the lift coefficients seen during translation are relatively constant. During rotation(Fig. 10E), the downstroke vortices are shed. After rotation, the wing moves back through its wake and new vortices are formed on the leading and trailing edges of the wing(Fig. 10F-I). These vortices remain attached to the wing during the upstroke and would be shed during the rotation at the beginning of the next stroke. In the case of larger insects(i.e. higher Re), lift is generated when the LEV remains attached and the trailing edge vortex is shed. When the trailing edge vortex remains attached, positive vorticity is not shed from the wing, and negative circulation around the wing is reduced (see the Discussion for an explanation). Finally, the strength of the wake is diminished compared with the larger Re case since viscous forces are relatively larger.

Another difference between the two cases is that vorticity dissipates relatively faster at lower Re. This can be seen by comparing the wake left by the downstroke in each case towards the end of the simulation. Any lift- or drag-altering effects produced when the wing moves through its wake will be diminished at lower Re. This wake capture effect should decrease gradually with decreasing Re.

Streamline plots for an Re of 16 were very similar to those for an Re of 8, and streamline plots for an Re of 64 were very similar to those for an Re of 128. This division would appear to be related to the transition seen behind steady plates and cylinders when the von Karman vortex street forms at an Re of ∼40. The simulation at an Re of 32 appears to be a borderline case. The streamline plots during the downstroke are very similar to those of an Re of 8. During the upstroke, the LEV begins to shed, and a von Karman vortex street might develop. Since the effective fluid velocity relative to the wing is larger during the upstroke (as the wing moves back through its wake), the effective Re would be transiently higher. This could account for some variation as the flow regime nears the transition Re.

To investigate the effects of Re on lift and drag generated at different angles of attack, we considered five angles at an Re of 8 and 128. In each case, the angle of attack during the downstroke was the same as the angle of attack on the upstroke. Changing the angle of attack also had the effect of changing the angle through which the wing was rotated and the angular velocity, since the duration of rotation was held constant for the different cases. All other kinematic parameters were the same as in the previous simulations.

Drag coefficients as a function of time (fraction of stroke) are plotted in Figs 11, 12. For both high and low Re, the drag coefficient increases with angle of attack. The drag coefficients are also substantially higher during periods of acceleration than during periods of constant translation in all cases. Drag coefficients reach their largest magnitudes shortly before and/or after changing sign during wing rotation. This effect is strongest at an angle of attack of 10°, because the wing rotates faster through larger angles than at higher angles of attack. In interpreting Figs 11 and 12 it is important to keep in mind that the pivot point is not in the center of the wing rotation but rather is located 0.2 chord lengths from the leading edge.

Lift coefficients as a function of time (fraction of stroke) are plotted in Figs 13, 14 for high and low Re, respectively. For both Re, the lift coefficients are greatest near an angle of attack of ∼40°. At Re=8,fluctuations in lift during translation are significantly lower than at Re=128. The lift coefficients are also larger when the wing accelerates than during translation in all cases. Lift coefficients reach their largest values during wing rotation. Similar to drag, this effect is strongest at an angle of attack of 10°. For both Re, lift drops significantly during the beginning of upstroke translation for low angles of attack. These lift coefficients approach downstroke values later during the upstroke.

To test for convergence, we ran two simulations: one at the mesh size used for all previous simulations and another at about half that mesh size. The first simulation used a 600×600 grid and the other used a 1200×1200 grid. Both simulations used the stroke kinematics described in Fig. 2 with a 45° angle of attack and an Re of 128. The resulting drag and lift coefficients as a function of dimensionless time are plotted in Figs 15, 16, respectively. The calculated lift and drag coefficients show good agreement, with small deviations during periods of wing acceleration and deceleration. The highest Re is shown for the convergence study because it is the most difficult case. Results at lower Re (not shown) yield better agreement.

In order to compare our simulation results with experimental data, we ran a simulation of a wing started almost impulsively from rest and translated at a constant speed over a distance of 7 chord lengths. The wing was accelerated from rest at a rate of 0.625 m s^-2 until it reached a translational speed of 0.10 m s^-1. This simulation modeled the experiments of Dickinson and Götz (1993)as closely as possible, using the same dimensions for the fluid domain (in this case, 1 m in length × 0.4 m in width) and the same chord length of the wing. Since this is a higher Re simulation than previous cases,the size of the fluid grid was increased to 1200×1200.

Figs 17, 18 compare the drag and lift coefficients at a 45° angle of attack with those measured by Dickinson and Götz (1993). In our simulations, lift oscillates with the alternate shedding of the leading and trailing edge vortices. Similar vortex shedding was observed by flow visualization in the Dickinson and Götz experiments. However, our oscillations in lift force are larger than those measured by Dickinson and Götz. Oscillations in drag coefficients in our simulations also correspond to the alternate shedding of the leading and trailing edge vortices but are twice the frequency of the oscillations in lift. This difference in the oscillation frequencies of the lift and drag forces is similar to what has been found for flow past cylinders in this Re range(Lai and Peskin, 2000). This difference can be explained by the fact that the shedding of either the leading or trailing edge vortices transiently reduces the drag force. However,the shedding of the LEV reduces lift while the shedding of the trailing edge vortex augments lift. The drag oscillations in our simulation are smaller in amplitude than the lift oscillations and match reasonably well with those measured by Dickinson and Götz.

The discrepancies between the results of our simulations and the experiments of Dickinson and Götz are unclear. Force oscillations in their experiment decrease significantly during the 7 chord translation. In our simulations, the amplitudes of force oscillations are relatively steady. Numerical error and experimental error probably account for some of the differences. Other differences might be explained in part by minor 3-D effects. While the Dickinson and Götz experiment was nearly 2-D, there were necessarily some edge effects at the span-wise ends of the wings. Other 3-D effects might include any span-wise flexing of the wing, although this effect would most likely be minor. Dickinson and Götz also found that the net force acting on the wing was approximately normal to the chord of the wing. This is not the case in our simulation, since oscillations in lift are larger than oscillations in drag, suggesting that viscous effects in our simulation are significant. This might not be entirely unreasonable. Vandenberghe et al. (2004)have shown that force can be generated tangent to the chord of a flat plate that oscillates in the direction normal to the chord of the plate, producing thrust. Moreover, alternate vortex shedding can generate large forces perpendicular to flow and tangent to the chord of a flat plate, causing`flutter' or auto-rotation in the direction tangent to the chord of the plate(Mittal et al., 2004; Skews, 1990).

Probably the most interesting result shown by these simulations is the aerodynamic transition observed between Re≥64 and Re≤32. The fundamental difference is that vortices are formed but not shed during translation at lower Re, but they are alternately shed during translation at higher Re. Similar transitions have been found in this intermediate range of Re for flow past a variety of shapes (Batchelor, 1967). This transition is significant because, below it, an important mechanism of lift generation may no longer apply. At some Re above 32, as the Re is increased, vortical asymmetry is produced when vortices are alternately shed. In 3-D insect flight, this asymmetry is manifested as an attached LEV throughout the translation of each half stroke and a shed trailing edge vortex. Such asymmetry leads to increased lift forces during translation (see discussion below). At some Re below 64, as the Re is reduced, vortical near-symmetry is produced when both the leading and trailing edge vortices remain attached to the wing during each half stroke in two dimensions. This symmetry would, in principle, reduce the amount of lift generated when compared with the asymmetric case in 3-D flight. Further research, however, is needed to see if this near-symmetry also occurs in three dimensions.

This aerodynamic transition between vortical symmetry and asymmetry is most likely related to similar transitions around the same Re seen in flow past cylinders and thrust generation in flapping flight. Between an Re of 4 and 40, the wake behind a cylinder consists of two symmetrical attached vortices, and no lift forces (or forces perpendicular to the flow) are produced. For Re above 40, vortices are alternately shed from each side of the cylinder, forming the well-known von Karman vortex street (Acheson, 1990; Batchelor, 1967). This vortical structure leads to alternating positive and negative forces on the cylinder perpendicular to the flow. Childress and Dudley(2004) describe a similar transition for thrust generation between an Re of 5 and 20. They considered the case of a wing flapping in a strictly vertical motion. Above some critical Re, this motion produces thrust (horizontal force). Vandenberhe et al. (2004) confirmed this transition in thrust production experimentally using an oscillating plate that was allowed to rotate perpendicular to the direction of the oscillations. These transitions in lift and thrust generation are most likely the result of a bifurcation in

A question that remains, however, is how well the 2-D models of insect flight at low Re apply to 3-D flight in tiny insects. There are several 3-D components that could be significant. First of all, actual wings are of finite span, whereas 2-D models assume infinite span. Secondly, the chord length of the wing varies with span, while the 2-D approximation assumes constant chord. Most importantly, the dorsal-ventral motion of the wings through translation is actually rotational (the wing is rotating at its root). At higher Re, these differences are significant. In 3-D flapping flight, alternate vortex shedding does not occur: the LEV remains attached to the wing until wing reversal and the trailing edge vortex is shed. This phenomenon appears to be robust for a range of Re leading to larger lift forces (Birch et al.,2004). Our 2-D simulations suggest that this vortical asymmetry would be lost at lower Re because the trailing edge vortex would not be shed. Such a loss of asymmetry in three dimensions would result in relatively lower lift forces for smaller insects. It is possible, however,that 3-D effects could induce the shedding of the trailing edge vortex for Re below the 2-D transition. Future work in three dimensions is necessary to verify this conclusion.

where F_L is the lift force on the body and F_D is the drag force on the body. These equations mean that lift and drag forces are proportional to the time rate of change of the total first moment of the vorticity field.

where|ω| is the absolute value of the vorticity. The magnitude of lift generated depends on the difference between time rate of change of the total first moment of positive vorticity and the time rate of change of the total first moment of negative vorticity. Due to the asymmetry in the vortical pattern behind the wing, positive vorticity is convected away from the wing at a greater rate than negative vorticity. Since total positive and negative vorticity in R_f must be equal (by the principle of vorticity conservation), this means that the total time rate of change of the first moment of positive vorticity is greater in magnitude than the total time rate of change of the first moment of negative vorticity. As a result,positive lift is produced.

For Re below 32, regions of negative and positive vorticity remain attached to the wing as the leading and trailing edge vortices (see Fig. 19B). The vortical pattern behind the wing is nearly symmetrical (true symmetry would occur at a 90° angle of attack). As a result, the difference between the time rate of change of the total first moment of vorticity in the leading and trailing edge vortices is reduced, and the lift coefficients are lower than for higher Re. At a 90° angle of attack, the time rate of change of the first moment of negative and positive vorticity would be balanced, producing no lift force.

List of symbols and abbreviations
 
• C_L

  lift coefficient

 
• C_D

  drag coefficient

 
• c_targ

  damping coefficient

 
• CFD

  computational fluid dynamics

 
• chord

  chord length of the wing

 
• DPIV

  digital particle image velocimetry

 
• f(r, t)

  force per unit area applied to the fluid by the wing

 
• f_beam(r, t)

  force per unit area applied to the fluid due to bending stiffness

 
• f_str(r, t)

  force per unit area applied to the fluid due to stretching stiffness

 
• f_targ(r, t)

  force per unit area applied by the target boundary

 
• F_D

  drag force per unit length

 
• F_L

  lift force per unit length

 
• F(x, t)

  force per unit volume acting on the fluid

 
• k_targ

  stiffness coefficient of the target boundary

 
• k_beam

  flexural stiffness coefficient of the wing

 
• k_str

  stiffness coefficient of the wing proportional to resistance to stretching

 
• l

  characteristic length of the wing

 
• LEV

  leading edge vortex

 
• M

  first moment of vorticity

 
• p(x, t)

  fluid pressure

 
• r

  Lagrangian position parameter

 
• R_∞

  two-dimensional infinite space

 
• R_f

  two-dimensional fluid space

 
• R_n

  region of negative vorticity

 
• R_o

  region of negligible vorticity

 
• R_p

  region of positive vorticity

 
• Re

  Reynolds number

 
• S

  surface area per unit length

 
• S

  two-dimensional region occupied by solid bodies

 
• U

  characteristic velocity

 
• u(x, t)

  fluid velocity

 
• V

  target translational velocity

 
• v(τ)

  translational velocity at dimensionless time τ

 
• t

  time

 
• x

  two-dimensional position vector

 
• X(r, t)

  Cartesian coordinate vector of the wing

 
• Y(r, t)

  Cartesian coordinate vector of the target boundary

 
• α

  angle of attack relative to the horizontal

 
• δ(x)

  two-dimensional delta function

 
• Δθ

  angular distance of rotation

 
• Δτtrans

  dimensionless duration of acceleration/deceleration

 
• Δτrot

  dimensionless duration of rotation

 
• θ

  angle between wing and positive x-axis

 
• θrot

  rotational constant

 
• \({\dot{{\theta}}}({\tau})\)

  angular velocity as a function of dimensionless time

 
• μ

  dynamic viscosity

 
• ρ

  fluid density

 
• τ

  dimensionless time

 
• τd

  dimensionless time of deceleration

 
• τfinal

  dimensionless duration of entire stroke

 
• τrot

  dimensionless time when rotation begins

 
• ν

  kinematic viscosity vorticity in a two-dimensional fluid

 
• |ω|

  absolute value of vorticity

 
• \({\bar{{\omega}}}\)

  flapping frequency

 
• ψ(x, t)

  stream function

We wish to thank Stephen Childress, Marvin Jones, Michael Dickinson and Will Dickson for many helpful conversations on insect flight and fluid dynamics. We also thank the reviewers for their insightful comments. This work was supported by the National Science Foundation under Grant Number DMS-9980069.