First published online July 26, 2004
Journal of Experimental Biology 207, 3073-3088 (2004)
Published by The Company of Biologists 2004
doi: 10.1242/jeb.01138
When vortices stick: an aerodynamic transition in tiny insect flight
Laura A. Miller* and
Charles S. Peskin
Courant Institute of Mathematical Sciences, New York University, 251
Mercer Street, New York, NY 10012, USA
*
Author for correspondence (e-mail:
millerl{at}cims.nyu.edu)
Accepted 10 June 2004
 |
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.
Key words: insect flight, Reynolds number, aerodynamics, computational fluid dynamics
|
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 mx1 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.
View larger version (19K):
[in this window]
[in a new window]
|
Fig. 1. A two-dimensional approximation of a three-dimensional stroke. The motion
of the wing is divided here into three stages: downstroke (A), rotation (B)
and upstroke (C). In reality, the rotational phase overlaps with the
downstroke and the upstroke. The wing moves approximately along a horizontal
plane. The center of rotation is 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) |
View larger version (14K):
[in this window]
[in a new window]
|
Fig. 2. Translational velocity and the angular velocity of the wing as a function
of dimensionless time for one stroke cycle. This motion was used for all
simulations unless otherwise stated. Note that the wing begins to rotate
during the first half stroke (or downstroke). Since most of the rotation
occurs at the end of the downstroke, this motion describes the case of
`advanced rotation'.
|
|
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.
View larger version (34K):
[in this window]
[in a new window]
|
Fig. 3. This diagram illustrates the numerical method used for these simulations:
the immersed boundary method. The fluid domain is represented as a Cartesian
grid. The boundary (wing) points are represented as red squares. These points
interact with the fluid and move at the local fluid velocity. The green
springs represent the bending and stretching stiffness of the boundary. The
desired motion of the wing is prescribed by the target points, which are shown
as blue circles. These points do not interact with the fluid and they move
according to the desired motion of the wing. They also apply a force to the
actual boundary via the target springs (shown in purple). The further
the actual boundary is from its target boundary, the larger the force applied
to the actual boundary.
|
|
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. ftarg(r, t) is
the force per unit area, ktarg is a stiffness coefficient,
ctarg 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. fbeam(r,
t) is the force per unit area, and kbeam is a
stiffness coefficient. Equation 14describes the force applied to the fluid as
a result of the resistance to stretching by boundary given as
fstr(r, t), where kstr 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 630x630. 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 CL is the lift coefficient,
CD is the drag coefficient, S is the surface area
per unit length of the model wing, U is the velocity of the boundary,
FD is the drag force per unit length,
FL 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
CD and CL 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.
View larger version (37K):
[in this window]
[in a new window]
|
Fig. 4. Drag coefficients are plotted as functions of time for one stroke cycle.
The arrows along the axis show the times at which streamline plots in Figs
9,
10 were drawn. The angles of
attack were chosen to produce a symmetric stroke. In all cases, the angle of
attack was 45°. Reynolds number (Re) was varied by changing the
translational velocity of the wing from 0.00375 to 0.06 m s-1. In
general, drag coefficients increase with decreasing Re. Maximum drag
forces occur during acceleration from rest at the beginning of the downstroke
and rotation at the end of the downstroke and at the beginning of the
upstroke.
|
|
View larger version (35K):
[in this window]
[in a new window]
|
Fig. 9. The plots show the streamlines of fluid flow around a flapping wing for one
stroke cycle starting from rest, at an Re of 128. The arrow on the
wing shows the direction of the normalized aerodynamic forces acting on the
wing. The angle of attack during pure translation was 45° for both the
downstroke and the upstroke. The maximum translational velocity was 0.06 m
s-1. The colors of the streamline reflect the value of the stream
function,  , along the streamlines. Red denotes the most positive values
and blue denotes the most negative values. During the downstroke, an attached
leading edge vortex (LEV) is initially formed while the trailing edge vortex
is shed (A,B). This corresponds to a growth in lift forces. In C, the LEV is
being shed while a new trailing edge vortex is formed. This corresponds to a
drop in lift. During rotation (D,E), the leading and trailing edge vortices
are shed. At the beginning of the upstroke (F), a new LEV is formed and a new
trailing edge vortex is formed and shed. This corresponds to an increase in
lift. In G, the LEV is shed and a new trailing edge vortex is formed. This
results in a drop in lift. Finally, a second leading edge vortex is formed and
the trailing edge vortex is shed, resulting in another lift peak (H,I).
|
|
View larger version (40K):
[in this window]
[in a new window]
|
Fig. 10. Streamline plots of fluid velocity around a flapping wing for one stroke
cycle starting from rest, at an Re of 8. The arrow on the wing shows
the direction of the normalized aerodynamic forces acting on the wing. The
angle of attack during pure translation was 45° for both the downstroke
and the upstroke. The maximum translational velocity was 0.00375 m
s-1. The colors of the streamline reflect the value of the stream
function, , along the streamlines. Red denotes the most positive values
and blue denotes the most negative values. Note that both the leading and the
trailing edge vortices remain attached to the wing except during stroke
reversal (A-D and F-I). This differs from the higher Re case, where
lift forces oscillate due to the alternate shedding of leading and trailing
edge vortices. Similar vortex dynamics were observed for Re up to
32.
|
|
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 Re32 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.
View larger version (8K):
[in this window]
[in a new window]
|
Fig. 8. The diagram above shows lift/drag averaged during steady downstroke
translation at a constant angle of attack of 45° plotted against
log10 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.
View larger version (32K):
[in this window]
[in a new window]
|
Fig. 11. Drag coefficients as functions of time are plotted for five angles of
attack for an Re of 128. For each simulation, the same angle of
attack was used on the downstroke and upstroke. The angles of attack for the
five simulations were 10°, 20°; 30°; 40° and 50°. The
maximum translational velocity of the wing was 0.06 m s-1. Drag
coefficients increase with increasing angle of attack. Maximum drag forces
occur during acceleration from rest at the beginning of the stroke and during
rotation at the end of the downstroke and at the beginning of the upstroke.
Drag forces are larger during rotation at lower angles of attack because the
distance over which the wing moves and its angular velocity are larger.
|
|
View larger version (34K):
[in this window]
[in a new window]
|
Fig. 12. Drag coefficients as functions of time are plotted for five angles of
attack for an Re of 8. For each simulation, the same angle of attack
was used on the downstroke and upstroke. The angles of attack for the five
simulations were 10°, 20°, 30°, 40° and 50°. The maximum
translational velocity of the wing was 0.00375 m s-1. Drag
coefficients increase with increasing angle of attack. Maximum drag forces
occur during acceleration from rest at the beginning of the stroke and during
rotation at the end of the downstroke and at the beginning of the upstroke.
Drag forces are larger during rotation at lower angles of attack because the
distance over which the wing moves and its angular velocity are larger.
|
|
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.
View larger version (32K):
[in this window]
[in a new window]
|
Fig. 13. Lift coefficients as functions of time are plotted for five angles of
attack for an Re of 128. For each simulation, the same angle of
attack was used on the downstroke and upstroke. The angles of attack for the
five simulations were 10°, 20°, 30°, 40° and 50°. The
maximum translational velocity of the wing was 0.06 m s-1. Lift
coefficients are greatest for an angle of attack near 40°. Maximum lift
forces occur during acceleration from rest at the beginning of each half
stroke and during rotation at the end of the downstroke and at the beginning
of the upstroke.
|
|
View larger version (29K):
[in this window]
[in a new window]
|
Fig. 14. Lift coefficients as functions of time are plotted for five angles of
attack for an Re of 8. For each simulation, the same angle of attack
was used on the downstroke and upstroke. The angles of attack for the five
simulations were 10°, 20°, 30°, 40° and 50°. The maximum
translational velocity of the wing was 0.00375 m s-1. During
translation, lift coefficients increase with increasing angle of attack in the
range of 10-40°, and lift coefficients for angles of attack of 40° and
50° are quite similar. Maximum lift forces occur during acceleration from
rest at the beginning of each half stroke and during rotation at the end of
the downstroke and beginning of the upstroke. Lift forces are larger during
rotation at lower angles of attack because the distance over which the wing
moves and its angular velocity are larger.
|
|
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 600x600 grid and the other used a
1200x1200 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.
View larger version (27K):
[in this window]
[in a new window]
|
Fig. 15. Drag coefficients for two mesh widths. The 600x600 grid size was used
for other simulations in this paper. Both simulations used the stroke
kinematics described in Fig. 2
at an Re of 128 and with an angle of attack of 45°.
|
|
View larger version (27K):
[in this window]
[in a new window]
|
Fig. 16. Lift coefficients for two mesh widths. The 600x600 mesh was used in
all other simulations in this paper. Both simulations used the stroke
kinematics described in Fig. 2
at an Re of 128 and with an angle of attack of 45°.
|
|
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 x 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 1200x1200.
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.
View larger version (27K):
[in this window]
[in a new window]
|
Fig. 17. Drag coefficients as a function of distance traveled for a wing started
from rest and translated seven chord lengths at a 45° angle of attack.
Dotted lines represent data collected during an experiment by Dickinson and
Götz (1993), and solid
lines represent the results of our two-dimensional simulation. Oscillations in
drag are smaller than those measured in lift and correspond to the alternate
shedding of the leading and trailing edge vortices.
|
|
View larger version (26K):
[in this window]
[in a new window]
|
Fig. 18. Lift coefficients as a function of distance traveled for a wing started
from rest and translated seven chord lengths at a 45° angle of attack.
Dotted lines represent data collected during an experiment of Dickinson and
Götz (1993), and solid
lines represent the results of our two-dimensional simulation. Lift
coefficients in both cases oscillate with the shedding of the leading and
trailing edge vortices. Note that lift forces in our simulation have stronger
oscillations than the forces 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 Re64 and
Re32. 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 Rf
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
[=(dv/dx)-(du/dy)], 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=[M1, M2]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 FL is the lift force on the body and
FD 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, Rf
(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 Rn be the
region of negative vorticity and Rp the region of
positive vorticity. Let Ro define the rest of
Rf with negligible vorticity. The total lift force
can then be calculated as follows:
 | (24) |
<
TOP
Summary
Introduction
