## SUMMARY

We have used computational fluid dynamics to study changes in lift
generation and vortex dynamics for Reynolds numbers (*Re*) between 8
and 128. The immersed boundary method was used to model a two-dimensional wing
through one stroke cycle. We calculated lift and drag coefficients as a
function of time and related changes in lift to the shedding or attachment of
the leading and trailing edge vortices.

We find that the fluid dynamics around the wing fall into two distinct
patterns. For *Re*≥64, leading and trailing edge vortices are
alternately shed behind the wing, forming the von Karman vortex street. For
*Re*≤32, the leading and trailing edge vortices remain attached to
the wing during each `half stroke'. In three-dimensional studies, large lift
forces are produced by `vortical asymmetry' when the leading edge vortex
remains attached to the wing for the duration of each half stroke and the
trailing edge vortex is shed. Our two-dimensional study suggests that this
asymmetry is lost for *Re* below some critical value (between 32 and
64), resulting in lower lift forces. We suggest that this transition in fluid
dynamics is significant for lift generation in tiny insects.

## Introduction

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.

## Materials and methods

Our 2-D numerical simulation of flight was constructed to be similar to the
physical experiment by Dickinson and Götz
(1993) but using a 2-D motion
similar to later experiments and numerical simulations
(Dickinson et al., 1999;
Sun and Tang, 2002). Dickinson
and Götz designed a single robotic wing to model flight similar to that
of *Drosophila melanogaster* and to understand better the aerodynamic
forces generated using flow visualization and direct force measurements. These
experiments were dynamically scaled such that the *Re* of the model was
approximately that of the flight of *D. melanogaster*. The *Re*
is a dimensionless variable that gives the ratio of inertial to viscous
forces:
1

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^{-1}
s^{-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.

The translational velocities throughout the stroke were constructed using a series of equations to describe each part of the stroke (see Fig. 2). The velocity during acceleration at the beginning of the downstroke is given by: 2 3

where τ is dimensionless time defined by equation 3, *t* is the
actual time, *chord* is the chord length of the wing, *V* is the
target translational velocity, *v*(τ) is the translational velocity
at dimensionless time τ, and Δτ_{trans} is dimensionless
duration of both the acceleration and deceleration phases of translation.
After acceleration, the wing travels with constant translational velocity
*V*. At the end of the downstroke, the deceleration of the wing is
given by:
4
5

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, τ_{d}
was 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 angle of attack relative to the horizontal axis was held constant
during the entire stroke except during the rotational phase at the end of the
downstroke and the beginning of the upstroke. Let α be the angle of
attack relative to the horizontal plane, and let θ be the angle between
the wing and the positive *x*-axis (the origin is defined as the
intersection of the wing with the *x*-axis). The angular velocity
during the rotational phase is given by:
6
7
where is the angular
velocity as a function of dimensionless time, θ_{rot} is a
constant determined by the distance of rotation and duration of the rotational
phase, Δτ_{rot} is the dimensionless duration of the
rotational phase, τ_{rot} is the dimensionless time when rotation
begins, and Δθ is the angular distance over which rotation occurs.
Unless otherwise noted, the value of θ during the following simulations
was 135° during the downstroke and 45° during the upstroke (note that
this corresponds to the same angle of attack, α=45°, in both cases).
Thus, Δθ was set to 90°. Also, Δτ_{rot} was
taken to be 3.48, and τ_{rot} to be 3. These values are similar to
those used by Dickinson et al.
(1999) and Sun and Tang
(2002) in the case of
`advanced rotation'.

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.

The equations of motion are as follows: 8 9

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.

The interaction equations between the fluid and the boundary are given by: 10 11

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*.

The immersed boundary equations are given by: 12 13 14 15

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 in the nonlinear term of the Navier-Stokes
equations was discretized as a skew symmetric operator to remove the effects
of numerical viscosity (Lai and Peskin,
2000). Essentially, the fluid equations are discretized on a
regular rectangular grid in the physical space of the position variable
**x**, and the boundary equations are discretized in a one-dimensional
space of the Lagrangian parameter *r*. The fluid domain is assumed to
be periodic. However, the periodicity in these simulations was broken by
including a stiff outer boundary near the edges of the domain. The dimensions
of the physical domain defined by the stiff outer boundary measure 1 m in
width and 1 m in length. The computational (periodic) domain was slightly
larger: 1.05 m in width and in length. The Eulerian fluid grid covering this
computational domain was 630×630. The immersed boundary (wing) was
discretized into 60 spatial steps. The stiffness coefficients were chosen to
reduce the deformation of the boundary to acceptable levels, and the damping
coefficient was chosen to provide light damping.

Lift and drag forces were calculated as a function of time by summing the forces that each immersed boundary point of the model wing applied to the fluid at each time step and taking the opposite sign of that value. Lift and drag coefficients were filtered to remove high frequency `noise' from the vibrations of the elastic boundary. This did not change the basic shape of the graphs. The lift and drag coefficients are defined as follows: 16 17

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, *S*
is 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.

## Results

### Changes in Re

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. 8
shows the average lift to drag ratio during translation for the downstroke as
a function of *Re*. Lift/drag increases with increasing
*Re*.

The aerodynamic basis of these *Re* changes may be seen by studying
Figs 9,
10. These figures are graphs
of the streamlines of the fluid flow around the wing for *Re* of 128
and 8 taken at 10 points during the stroke. These points in time are shown by
arrows in Figs 4,
6. The streamlines are curves
that have the same direction as the instantaneous fluid velocity,
**u**(**x**, *t*), at each point. They were drawn by making a
contour map of the stream function, since the stream function is constant
along streamlines. The stream function ψ(**x**, *t*) in 2-D is
defined by the following equations:
18

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-C
shows 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*.

### Changes in angle of attack

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.

### Convergence test

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.

### Comparison with experimental data

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).

## Discussion

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
,
where ρ is the density of the fluid,
is the flapping frequency,
*L* is the body length, and μ is the viscosity of the fluid
(Childress and Dudley,
2004).

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.

To understand how vortical asymmetry leads to lift generation, we first
present the general aerodynamic theory for viscous flows given by Wu
(1981). Consider a 2-D viscous
fluid initially at rest with an immersed solid body also initially at rest in
an infinite space, **R**_{∞}. Let **R**_{f}
define the space occupied by the fluid, and **S** define the space occupied
by the solid. Since the total vorticity in **R**_{∞} is
initially zero, by the principle of total vorticity conservation the total
vorticity in **R**_{∞} is zero for all time:
19

where **x** is the position vector x=(*x, y*)^{T},ω
is the vorticity in a two-dimensional flow
[ω=(d*v*/d*x*)-(d*u*/d*y*)], and the fluid
velocity is given as **u**(**x**, *t*)=[*u*(**x**,
*t*), *v*(**x**, *t*)]^{T}. Note that this
principle is only true because we are considering vorticity in the total space
occupied by both the solid and the fluid. Wu
(1981) showed that the
aerodynamic force exerted on the solid body is given as follows:
20
21

where **M**=[*M*_{1}, *M*_{2}]^{T}
is the first moment of the vorticity field, ρ is the density of the fluid
and **S** is the region occupied by the solid. The second term in equation
20 is an inertial term for the body. During periods of constant translation,
this term goes to zero and we have the following equations:
22
23

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.

Consider the case for a wing translated from rest with an attached LEV and
a shed trailing edge vortex in a region of fluid, **R**_{f}
(as shown in Fig. 19A). In the
following discussion, the coordinate axis moves with the center point of the
boundary (so that in this frame of reference the boundary is at rest and the
fluid moves past it). Positive flow moves from the left to the right. An
attached region of negative (clockwise) vorticity forms along the leading edge
of the wing. Positive (counterclockwise) vorticity is shed from the wing in
the form of a starting vortex and wake. Let **R**_{n} be the
region of negative vorticity and **R**_{p} the region of
positive vorticity. Let **R**_{o} define the rest of
**R**_{f} with negligible vorticity. The total lift force
can then be calculated as follows:
24

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
- 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
- flapping frequency
- ψ(
**x**,*t*) - stream function

## ACKNOWLEDGEMENTS

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.

- © The Company of Biologists Limited 2004