When vortices stick: an aerodynamic transition in tiny insect flight
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 find that the fluid dynamics around the wing fall into two distinct
patterns. For Re64, 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 |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
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 |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
![]() | (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.
|
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) |
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. 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 |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
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 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.
|
|
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.
|
|
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 |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
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) |
|
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 Rf 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 |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
![]() |
Acknowledgments |
---|
![]() |
References |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
Acheson, D. J. (1990). Elementary Fluid Dynamics. Oxford: Clarendon Press.
Batchelor, G. K. (1967). An Introduction to Fluid Dynamics. Cambridge: Cambridge University Press.
Birch, J. M. and Dickinson, M. H. (2003). The
influence of wing-wake interactions on the production of aerodynamic forces in
flapping flight. J. Exp. Biol.
206,2257
-2272.
Birch, J. M., Dickson, W. B. and Dickinson, M. H.
(2004). Force production and flow structure of the leading edge
vortex on flapping wings at high and low Reynolds numbers. J. Exp.
Biol. 207,1063
-1072.
Childress, S. and Dudley, R. (2004). Transition from ciliary to flapping mode in a swimming mollusc: flapping flight as a bifurcation in Re_omega. J. Fluid. Mech. 498,257 -288.[CrossRef]
Cloupeau, M., Devillers, J. F. and Devezeaux, D. (1979). Direct measurements of instantaneous lift in desert locust: comparison with Jensen's experiments on detached wings. J. Exp. Biol. 80,1 -15.
Dickinson, M. H. (1994). The effects of wing
rotation on unsteady aerodynamic performance at low Reynolds numbers.
J. Exp. Biol. 192,179
-206.
Dickinson, M. H. and Götz, K. G. (1993).
Unsteady aerodynamic performance of model wings at low Reynolds numbers.
J. Exp. Biol. 174,45
-64.
Dickinson, M. H., Lehmann, F.-O. and Sane, S. P.
(1999). Wing rotation and the aerodynamic basis of insect flight.
Science 284,1954
-1960.
Dudley, R. (2000). The Biomechanics of Insect Flight: Form, Function, Evolution. Princeton: Princeton University Press.
Ellington, C. P. (1984a). The aerodynamics of hovering insect flight. I. The quasi-steady analysis. Phil. Trans. R. Soc. Lond. B 305,1 -15.
Ellington, C. P. (1984b). The aerodynamics of hovering insect flight. III. Kinematics. Phil. Trans. R. Soc. Lond. B 305,41 -78.
Ellington, C. P. (1995). Unsteady aerodynamics of insect flight. In Biological Fluid Dynamics, vol.49 (ed. C. P. Ellington and T. J. Pedley), pp.109 -129. Cambridge: The Company of Biologists Ltd.
Ellington, C. P. (1999). The novel aerodynamics
of insect flight: applications to micro-air vehicles. J. Exp.
Biol. 202,3439
-3448.
Ellington, C. P., van den Berg, C., Willmont, A. P. and Thomas, A. L. R. (1996). Leading edge vortices in insect flight. Nature 348,626 -630.[CrossRef]
Horridge, G. A. (1956). The flight of very small insects. Nature 178,1334 -1335.
Lai, M.-C. and Peskin, C. S. (2000). An immersed boundary method with formal second order accuracy and reduced numerical viscosity. J. Comp. Phys. 160,705 -719.[CrossRef]
Landau, L. D. and Lifshitz, E. M. (1959). Fluid Mechanics. London: Pergamon Press.
Lighthill, M. J. (1973). On the Weis-Fogh mechanism of lift generation. J. Fluid Mech. 60, 1-17.
Lighthill, M. J. (1975). Mathematical Biofluiddynamics. Philadelphia: SIAM.
Liu, H., Ellington, C. P., Kawachi, C., van den Berg, C. and
Willmott, A. P. (1998). A computational fluid dynamic study
of hawkmoth hovering. J. Exp. Biol.
201,461
-477.
Mittal, R., Veeraraghavan, S. and Udaykumar, H. S. (2004). Flutter, tumble, and vortex induced autorotation. Theor. Comp. Fluid Dynamics 17,165 -170.[CrossRef]
Osborne, M. F. M. (1951). Aerodynamics of flapping flight with application to insects. J. Exp. Biol. 28,221 -245.[Medline]
Peskin, C. S. (2002). The immersed boundary method. Acta Numerica 11,479 -517.[CrossRef]
Peskin, C. S. and McQueen, D. M. (1996). Fluid dynamics of the heart and its valves. In Case Studies in Mathematical Modeling - Ecology, Physiology, and Cell Biology (ed. H. G. Othmer, F. R. Adler, M. A. Lewis and J. C. Dallon), pp.309 -337. New Jersey: Prentice Hall.
Ramamurti, R. and Sandberg, W. C. (2002). A
three-dimensional computational study of the aerodynamic mechanism of insect
flight. J. Exp. Biol.
205,2997
-3008.
Sane, S. P. and Dickinson, M. H. (2002). The
aerodynamic effects of wing rotation and a revised quasi-steady model of
flapping flight. J. Exp. Biol.
205,1087
-1096.
Skews, B. W. (1990). Autorotation of rectangular plates. J. Fluid. Mech. 165,247 -272.
Sotavalta, O. (1947). The flight-tone (wing stroke frequency) of insects. Acta Entomol. Fenn. 4, 1-117.
Spedding, G. R. and Maxworthy, T. (1986). The generation of circulation and lift in a rigid two-dimensional fling. J. Fluid. Mech. 165,247 -272.
Sun, M. and Tang, J. (2002). Unsteady
aerodynamic force generation by a model fruit fly wing in flapping motion.
J. Exp. Biol. 205,55
-70.
Thom, A. and Swart, P. (1940). Forces on an airfoil at very low speeds. J. Roy. Aeronaut. Soc. 44,761 -770.
Thompson, D. A. W. (1917). On Growth and Form. Cambridge: Cambridge University Press.
Usherwood, J. R. and Ellington, C. P. (2002).
The aerodynamics of revolving wings I. Model hawkmoth wings. J.
Exp. Biol. 205,1547
-1564.
van den Berg, C. and Ellington, C. P. (1997a). The three-dimensional leading-edge vortex of a `hovering' model hawkmoth. Phil. Trans. R. Soc. Lond. B 352.
van den Berg, C. and Ellington, C. P. (1997b). The vortex wake of a `hovering' model hawkmoth. Phil. Trans. R. Soc. Lond. B 352,317 -328.[CrossRef]
Vandenberghe, N., Zhang, J. and Childress, S. (2004). Symmetry breaking leads to forward flapping flight. J. Fluid. Mech. 506,147 -155.[CrossRef]
Walker, J. A. (2002). Rotational lift: something different or more of the same? J. Exp. Biol. 205,3783 -3792.[Medline]
Wang, Z. J. (2000a). Shedding and frequency selection in flapping flight. J. Fluid. Mech. 410,323 -341.[CrossRef]
Wang, Z. J. (2000b). Two dimensional mechanism for insect hovering. Phys. Rev. Lett. 85,2216 -2218.[CrossRef][Medline]
Weis-Fogh, T. (1973). Quick estimates of flight fitness in hovering animals, including novel mechanisms for lift production. J. Exp. Biol. 59,169 -230.
Weis-Fogh, T. (1975). Unusual mechanisms for the generation of lift in flying animals. Sci. Am. 233, 80-87.
Willmott, A. P., Ellington, C. P. and Thomas, A. L. R. (1997). Flow visualization and unsteady aerodynamics in the flight of the hawkmoth Manduca sexta. Phil. Trans. R. Soc. Lond. B 352,303 -316.[CrossRef]
Wu, J. C. (1981). Theory for aerodynamic force and moment in viscous flows. AIAA J. 19,432 -441.
Wu, J. H. and Sun, M. (2004). Unsteady
aerodynamic forces of a flapping wing. J. Exp. Biol.
207,1137
-1150.
Zanker, J. M. and Götz, K. G. (1990). The wing beat of Drosophila melanogaster. II. Dynamics. Phil. Trans. R. Soc. Lond. B 327,19 -44.