A three-dimensional computational study of the aerodynamic mechanisms of insect flight
Ravi Ramamurti* and
William C. Sandberg
Laboratory for Computational Physics and Fluid Dynamics, Naval
Research Laboratory, Washington, DC 20375, USA
*
e-mail:
ravi{at}lcp.nrl.navy.mil
Accepted 4 March 2002
 |
Summary
|
---|
A finite element flow solver was employed to compute unsteady flow past a
three-dimensional Drosophila wing undergoing flapping motion. The
computed thrust and drag forces agreed well with results from a previous
experimental study. A grid-refinement study was performed to validate the
computational results, and a grid-independent solution was achieved. The
effect of phasing between the translational and rotational motions was studied
by varying the rotational motion prior to the stroke reversal. It was observed
that, when the wing rotation is advanced with respect to the stroke reversal,
the peak in the thrust forces is higher than when the wing rotation is in
phase with the stroke reversal and that the peak thrust is reduced further
when the wing rotation is delayed. As suggested by previous authors, we
observe that the rotational mechanism is important and that the combined
translational and rotational mechanisms are necessary to describe accurately
the force time histories and unsteady aerodynamics of flapping wings.
Key words: unsteady mechanism, incompressible fluid, flapping foil, unstructured mesh, insect, flight, Drosophila
 |
Introduction
|
---|
Flapping foil propulsion has received considerable attention in the past
few years as an alternative to the propeller. This mode of propulsion, which
involves no body undulation, has many applications, such as propulsion of
submersibles, maneuvering and flow control, that are of interest to the
hydrodynamic community and unconventional aerodynamics of micro aerial
vehicles (MAVs) and the study of aircraft flutter that are of interest to the
aerodynamic community.
Flapping foil propulsion is also important in the area of biofluid dynamics
for the study of propulsion in insects, birds and certain aquatic animals.
Flying animals generate lift and thrust as a consequence of the interaction
between the flapping motions of the wings and the surrounding air. These
animals also perform maneuvers involving rapid plunging and pitching motions.
Conventional steady-state theories do not predict sufficient forces to meet
those required for flight (Ellington,
1984
). Therefore, we need to understand the unsteady aerodynamics
of flapping wings undergoing highly three-dimensional motions with widely
varying geometries.
Experimental work on two-dimensional flapping foils has been carried out by
Anderson (1996
) and Freymuth
(1999
). Computational studies
have been performed by Jones and Platzer
(1997
) and Ramamurti and
Sandberg (2001
). While
two-dimensional wing section investigations can yield useful insights on the
coupled pitching and heaving dynamics, nothing can be learned concerning the
influence of spanwise flow. It is therefore essential to carry out
computations for actual three-dimensional insect wings. Ramamurti et al.
(1996
) computed
three-dimensional unsteady flow past moving and deforming geometries in a
simulation of a swimming tuna with caudal fin oscillation.
The three-dimensional wing strokes of insects can be divided into two
translational phases and two rotational phases. During the translational
phases, the upstroke and the downstroke, the wings move through the air with
high angles of attack, and during the rotational phases, pronation and
supination, the wings rotate rapidly and reverse direction. Dickinson et al.
(1999
) studied the effects of
translational and rotational mechanisms of the wing in Drosophila
melanogaster. They directly measured the forces produced by flapping
wings and explained the aerodynamics of insect flight by interactions between
three unsteady flow mechanisms. The `delayed' stall mechanism is a
translational mechanism which, in two dimensions, produces high lift in the
initial phases of translation until eventual flow separation; in three
dimensions, the spanwise flow effectively prevents stall. Rotational
circulation and wake capture are rotational mechanisms that depend mainly on
the pronation and supination of the wing during stroke reversal. Walker and
Westneat (1997
) have studied
experimentally the kinematics of fin motion in a fish, the bird wrasse
Gomphosus varius. Liu and Kawachi
(1998
) numerically investigated
the flow over a hovering hawkmoth, Manduca sexta. They reported the
presence of a spiral leadingedge vortex to which they attributed the lift
enhancement mechanism. They validated their results by comparing the
computational streak lines found in a two-dimensional hovering airfoil with
the experimental smoke visualization, but they did not directly compare
instantaneous forces.
Here, we extend the two-dimensional pitching and heaving airfoil
computations to three dimensions. This study will address the role of the
rotational motion in detail. Also, the role of the leading edge vortex and the
interaction between the axial flow and this leading edge vortex are
investigated. The primary objectives were (i) to validate the
three-dimensional computations by comparing the forces with the experimental
results of Dickinson et al.
(1999
) and performing
gridrefinement studies, to verify the hypothesis of Dickinson et al.
(1999
) that rotational
mechanisms of the wing form the basis by which the insect modulates the
magnitude and direction of the forces during flight and (ii) to provide data
on the forces, moments and power required for the development of a robotic
fly.
To this end, computations are performed for various phase angles between
the rotation and translation motions, and the time history of the unsteady
forces is compared with the experimental results. The flow solver we employ is
a finiteelement-based incompressible flow solver based on simple, low-order
elements. The simple elements enable the flow solver to be as fast as
possible, reducing the overheads in building element matrices, residual
vectors, etc. The governing equations are written in Arbitrary Lagrangian
Eulerian form, which enables flow with moving bodies to be simulated. The
details of the flow solver, the rigid body motion and adaptive remeshing are
given by Ramamurti et al.
(1995
) and are summarized
below.
 |
Materials and methods
|
---|
The incompressible flow solver
The governing equations employed are the incompressible NavierStokes
equations in Arbitrary Lagrangian Eulerian (ALE) formulation. They are written
as:
 | (1) |
 | (2) |
 | (3) |
Here, p denotes pressure, and va=v-w, the
advective velocity vector, where v is flow velocity and w is
mesh velocity and the material derivative is with respect to the mesh velocity
w. Both the pressure p and the stress tensor
have been
normalized by the (constant) density
and are discretized in time
t using an implicit time-stepping procedure. It is important for the
flow solver to be able to capture the unsteadiness of a flow field. The
present flow solver is built as time-accurate from the onset, allowing local
time stepping as an option. The resulting expressions are subsequently
discretized in space using a Galerkin procedure with linear tetrahedral
elements. To be as fast as possible, the overheads in building element
matrices, residual vectors, etc., should be kept to a minimum. This
requirement is met by employing simple, low-order elements that have all the
variables (v,p) at the same location. The resulting matrix
systems are solved iteratively using a preconditioned conjugate gradient
algorithm (PCG), as described by Martin and Löhner
(1992
). The flow solver has
been successfully evaluated for both two-dimensional and three-dimensional
laminar and turbulent flow problems by Ramamurti and Löhner
(1992
) and Ramamurti et al.
(1994
).
To carry out computations of the flow about oscillating and deforming
geometries, one needs to describe grid motion on a moving surface, i.e. to
couple the moving surface grid to the volume grid. The volume grid in the
proximity of the moving surface is then remeshed to eliminate badly distorted
elements. The velocity of the mesh is obtained in a manner so as to reduce
this distortion. A detailed description of the various mesh movement
algorithms is given in Ramamurti et al.
(1994
). In that study,
smoothing of the coordinates was employed for the mesh movement with a
specified number of layers of elements that move rigidly with the wing. In
two-dimensional studies (Ramamurti and
Sandberg, 2001
), the grid showed that the elements at the edge of
the rigid layers were very distorted after one cycle of oscillation. This is
due to a residual mesh velocity that is present as a result of the
non-convergence of the mesh velocity field. This will appear whether a
spring-analogy or a Laplacian-based smoothing is used.
To reduce the distortion of the mesh, the coordinates at the new time
xn+1 were obtained as a weighted average of the original
grid point location at time t=0(x0) and the location of
the point as if it moved rigidly with the body
(xrigidn+1):
 | (4) |
where the weighting function (f) is a simple linear function based on
the distance from the center of rotation r, and is given by:
 | (5) |
 | (6) |
and
 | (7) |
The mesh velocity is then obtained using:
 | (8) |
The values of rmin and rmax used in
this study are 2.0 and 10.0, respectively. To reduce the computational effort,
the experimental arrangement of Dickinson et al.
(1999
) was approximated by
introducing a symmetry plane. Because of the proximity of the wing at the
beginning of the downstroke and the rotation of the wing during the pronation
phase, the normal component mesh velocity of the points on the symmetry plane
can become non-zero. This would result in the points being pulled away from
the symmetry plane. To avoid this problem, the points on the symmetry plane
are allowed to glide along this plane. A similar technique has been employed
for the simulation of torpedo launch from a submarine by Ramamurti et al.
(1995
) where the gap between
the launch tube and the torpedo was small.
 |
Results and discussion
|
---|
The configuration of the hovering Drosophila melanogaster is shown
in Fig. 1A. The coordinate
system (x,y,z) is fixed to the body with the x coordinate
normal to the stroke plane. During the translational phases (upstroke and
downstroke), the wing moves from close to the y axis through an angle
, the wingbeat amplitude. The flapping wing configuration used in the
flow simulations is shown in Fig.
1B. This is based on the experimental arrangement of Dickinson et
al. (1999
).
Fig. 1B shows the position of
the wing at three different times during the flapping cycle. The coordinate
system (x',y',z') is fixed to the wing, and the
wing rotates about the z' axis throughout the cycle. The
planform of the wing is derived from the Drosophila wing and is 25 cm
long and 3.2 mm thick. The experimental apparatus consisted of two wings
immersed in a tank of mineral oil. The viscosity of the oil, the length of the
wing and the frequency of the flapping motion were chosen to match the
Reynolds number (Re) of a typical Drosophila melanogaster,
approximately 136. The Re for the present calculations is defined on
the basis of the mean chord of the wing
(6.7 cm) and the mean wing-tip
velocity Ut (ignoring the forward velocity), as follows:
 | (9) |
where
,
Ut=2
nR, R is the wing length (25 cm),
AR is the aspect ratio of the wing, n is the frequency of
flapping motion,
is the wingbeat amplitude (peak to peak, in rad) and v
is the kinematic viscosity (115
cSt=115x10-6m2s-1).

View larger version (16K):
[in this window]
[in a new window]
|
Fig. 1. (A) Schematic diagram of a hovering Drosophila showing the
orientation of the x,y,z coordinate system. (B) Schematic diagram of
the flapping Drosophila wing. The position of the wing is shown at
three different times during the flapping cycle. The coordinate system
(x',y',z') is fixed to the wing, and the wing
rotates about the z' axis throughout the cycle. R,
wing length; , wingbeat amplitude.
|
|
The kinematics of the wing motion is obtained from the experiments of
Dickinson et al. (1999
). The
angles of rotation about the x axis (roll) and the z'
axis (pitch) are shown in Fig.
2A. The motion of the wing is prescribed using these two angles.
Fig. 2B shows the translational
velocity of the wing tip and the rotational (angular) velocity of the wing.
Three different phases between the translational and rotational motions were
used. In the `advanced' case, wing rotation precedes stroke reversal by 8 % of
the wingbeat cycle; `symmetrical' wing motion is where the wing rotation
occurs symmetrically with respect to stroke reversal; in `delayed' wing
motion, rotation is delayed by 8 % with respect to stroke reversal. Wingbeat
amplitude is 160°, flapping frequency is 0.145 Hz and the angle of attack
at midstroke is approximately 40° during both upstroke and downstroke.

View larger version (24K):
[in this window]
[in a new window]
|
Fig. 2 . Kinematics of the flapping wing. (A) Angle of rotation of the wing about
the x axis (roll), and the z' axis (pitch) for three
different phases between wing rotation and stroke reversal. (B) Translational
velocity of the wing tip and rotational (angular) velocity of the wing for
three different phases.
|
|
Symmetrical case
The flow solver described here is employed to compute the flow past the
Drosophila wing undergoing translation and rotation. First, an
inviscid solution was obtained using a grid consisting of 178 219 points and
965 877 tetrahedral elements. An initial steady-state solution was obtained in
1500 time steps. The unsteady solution using the prescribed kinematics
(Fig. 2) is then obtained. The
surface pressure on the wing is integrated to obtain the forces on the wing
along the three axes (Fx, Fy,
Fz). The thrust T and the drag D forces
are then computed as T=-Fx and
D=
(Fy2+Fz2),
respectively. These forces are compared with those obtained from the
experiments of Dickinson et al.
(1999
).
The unsteady computation was carried out for five cycles of oscillation.
Fig. 3 shows the thrust and
drag forces during one cycle of the wingbeat and compares the values with
those of Dickinson et al.
(1999
). The present
computations capture the peak forces well. The mean thrust force is
approximately 0.318 N, and the mean thrust coefficient
is 1.317. The
mean drag force is 0.375 N and the drag coefficient
is 1.55. The
force coefficients were obtained using the following non-dimensionalization:
 | (10) |
 | (11) |
where
and
are the mean thrust and drag forces,
respectively, and r22(s) is the second
moment of the dimensionless area of the wing (0.40). The variation of these
forces during the translational phase of the wing is also predicted correctly,
but the magnitude of the thrust force during the downstroke is higher than
that of Dickinson et al.
(1999
). The kinematics is
symmetrical between the up- and downstroke, so the resultant force should also
be symmetrical. That this is not the case may be due the mechanical play in
the experimental arrangement, as suggested by M. H. Dickinson (personal
communication). To understand the different mechanisms occurring during the
wingbeat cycle, we can divide the cycle into two rotational and two
translational phases. The rotational phase near the beginning of the
downstroke (pronation) occurs between time t0 and t3
(Fig. 3A). Thrust decreases
between t0 and t1, then increases until t2. This
behavior can be explained by a rotational mechanism. The wing continuously
rotates in a counterclockwise direction producing a circulation pointing
nearly along the +y direction. Between t0 and t1,
the wing is translating in the -z direction, resulting in a force
pointing in the -x direction, thus producing a peak in thrust at
t0. If a rotational mechanism alone were present, the thrust should
continue to decrease until t3; in fact, the thrust force increases
between t1 and t2. This happens after the wing changes
direction at the start of each half-stroke. Dickinson et al.
(1999
) attribute this increase
in thrust to a wake-capture mechanism in which the wing passes through the
shed vorticity of the previous stroke.

View larger version (23K):
[in this window]
[in a new window]
|
Fig. 3 . Time history of thrust (A) and drag (B) forces during one wingbeat. The
red lines are from the present study; the blue lines are from Dickinson et al.
(1999 ). The numbered broken
lines in A refer to times discussed in the text.
|
|
The position of the wing at the beginning of the downstroke is shown in
Fig. 4A. The chord in this case
of symmetrical rotation is aligned with the x axis at this instant.
We found a separation bubble attached to the leading edge during the interval
t0-t1. Particles were released from a rake of rectangular
grid of points in a plane 0.8 mm away from the bottom surface of the wing.
Using the instantaneous velocity field, the positions of these particles were
obtained by integrating the velocity at these rake points in both the positive
and negative velocity directions until the length of the traces exceeded a
specified length or the particles ended on a solid boundary or exited the
computational domain. These instantaneous particle traces are colored
according to the magnitude of velocity at that location.
Fig. 5 shows the leading edge
vortex with vorticity oriented in the +y direction. This leading edge
vortex was created at the end of the preceding upstroke. This vortex is
located below the wing and is rotating in the counterclockwise direction,
which can be seen from the velocity vectors shown in
Fig. 6A. A possible explanation
for the increase in thrust between t1 and t2 is that the
wing moving through this wake benefits from the shed vorticity. As the wing
moves through this vortex during the downstroke, it produces a stagnation
region at the bottom of the wing near the z' axis
(Fig. 5), resulting in an
increase in thrust.

View larger version (14K):
[in this window]
[in a new window]
|
Fig. 4 . (A) Position of the wing (red outline) at the beginning of the
downstroke. The wing chord is aligned with the x axis of the
x,y,z coordinate system at this instant. The orientation of the
y=10 cm plane for which the velocity vectors are shown in
Fig. 6 is indicated. (B)
Position of the wing at t=12.5 s. The orientation of the two planes
for which velocity vectors are shown in
Fig. 7 is indicated.
|
|

View larger version (74K):
[in this window]
[in a new window]
|
Fig. 5 . Instantaneous particle traces at the beginning of the downstroke.
Particles were released from a rake of rectangular grid of points in a plane
0.8 mm away from the bottom surface of the wing. Using the instantaneous
velocity field, the positions of these particles were obtained by integrating
the velocity at these rake points until the length of the traces exceeded a
specified length or the particles ended on a solid boundary or exited the
computational domain. These particle traces are colored according to the
magnitude of the velocity (in cm s-1) at that location. A leading
edge vortex is seen rotating in the counterclockwise direction, and a
stagnation line is shown near the z' axis of rotation (dark
blue traces).
|
|

View larger version (37K):
[in this window]
[in a new window]
|
Fig. 6. Velocity vectors near the leading edge at three instants at the beginning
of the downstroke on a plane y=10 cm (see
Fig. 4A for the orientation of
the plane relative to wing). L.E., leading edge. The vectors are colored
according to the magnitude of absolute velocity (cm s-1) and are of
constant length. The flow separates at the leading edge (A) and the separation
point moves down (B,C) as the downstroke is continued.
|
|

View larger version (58K):
[in this window]
[in a new window]
|
Fig. 7. Velocity vectors (in cm s-1; see color scale) during the
translation phase of the downstroke (between t3 and t4 in
Fig. 3A). (A,B,C) Velocity
vectors near the trailing edge on the xz plane at y=10 cm;
(D,E,F) velocity vectors near the wing tip on the yz plane at
x=0 cm (see Fig. 4B
for the orientation of planes relative to the wing). T.E., trailing edge.
|
|
At the beginning of the downstroke (t1), the flow separates at the
leading edge and reattaches on the bottom of the wing as shown in
Fig. 6A. As the wing continues
to move down, between t1 and t2, the separation point of
this bubble moves back along the wing chord
(Fig. 6). During the interval
t2-t3, we observe a trailing edge separation bubble forming
(see Fig. 7A-C). A similar
separation region forms at the wing tip, as can be seen in
Fig. 7D-F.
Fig. 4B shows the position of
the wing at t=12.5 s. Thrust production reaches a local minimum
around this instant (Fig. 3).
In Fig. 7A,D, a large
recirculation region can be seen in the wake of the wing. This separated flow
from the wing tip and trailing edge will result in a higher pressure on the
upper surface of the wing and, hence, a reduction in the thrust. During the
interval t1-t3, the magnitude of the translational
acceleration of the wing decreases while that of the angular acceleration
increases (symmetrical phase, Fig.
8).

View larger version (38K):
[in this window]
[in a new window]
|
Fig. 8. Translational and angular accelerations of the wing for three different
phases between wing rotation and stroke reversal. The broken vertical lines
refer to the sections of the wing stroke identified in
Fig. 3A.
|
|
Between t2 and t3, the magnitude of the translational
acceleration is large enough to overcome the rotational effect and, when the
angular acceleration becomes large enough, the rotational mechanism dominates,
resulting in the observed reduction in thrust until t3. Between
t3 and t4, the translational effect should result in a
constant thrust because the translational acceleration is almost constant
during this period. The rotational effect produces an increase in thrust
between t3 and t4, with a plateau in the middle, which
occurs when the trailing edge vortex is shed. Similar trends are observed
during the supination phase prior to the beginning of the upstroke
(t4-t5) and at the beginning of the upstroke
(t5-t7).
Fig. 9 shows the
instantaneous traces or streaklines of particles released 3.0 cm above or 3.5
cm below the wing in a plane parallel to the wing and near the leading edge.
In Fig. 9A, a wing tip vortex
can be seen, but no leading edge vortex is visible above the wing surface. A
leading edge vortex spinning in the counterclockwise direction is found below
the wing surface (Fig. 9B).
Particle traces near the mid stroke are shown in
Fig. 10A. At this instant, the
wing rotation axis z' is aligned with the body coordinate
z (see Fig. 1B). Here,
we can see the beginnings of the leading edge vortex on the upper surface of
the wing that is also shown by the velocity vectors in the xy plane
at z=20 cm (Fig.
10B). Fig. 10C
shows the pressure contours on the upper surface of the wing. The pressure is
non-dimensionalized with respect to the dynamic head,

Ut2, where
is the density of
the mineral oil (880 kg m-3) used in the experiments of Dickinson
et al. (1999
). A region of
constant pressure is present from the root of the wing up to approximately 60%
of the span and extending to the trailing edge.
Fig. 11 shows the particle
traces prior to the end of the downstroke. Here, the leading edge vortex is
clearly seen on the upper surface of the wing.

View larger version (32K):
[in this window]
[in a new window]
|
Fig. 9. Particle streaklines produced during the downstroke at t=12.93s.
The particles were release from a rake of rectangular grid on a plane 3.0 cm
above (A) and 3.5 cm below (B) the wing parallel to the wing and near the
leading edge. The absolute velocities of the particles (in cm s-1)
are given in the color scale.
|
|

View larger version (39K):
[in this window]
[in a new window]
|
Fig. 10. Flow patterns during the middle of the downstroke at t=13.79 s.
(A) Particle traces (streaklines). The particles were released on a plane
parallel to the wing and 3.9 cm below it, and velocity vectors are colored
according to the magnitude of absolute velocity (in cm s-1). (B)
Velocity vectors on the xy plane at z=20 cm. Velocity
vectors are colored according to the magnitude of absolute velocity (in cm
s-1) and are of constant length. (C) Pressure contours. Pressure is
non-dimensionalized with respect to the dynamic head.
|
|

View larger version (38K):
[in this window]
[in a new window]
|
Fig. 11. Particle traces prior to the end of the downstroke, t=14.65 s. The
particles were release from a rake of rectangular grid on a plane 3.2 cm above
the wing parallel to the wing and near the leading edge. Traces are colored
according to the magnitude of absolute velocity (in cm s-1).
|
|
Grid-refinement study
To assess the sensitivity of our computational results, we carried out a
grid-refinement study. The resolution of the grid in the vicinity of the wing
was doubled. The computations were carried out using a grid consisting of
approximately 238x103 points and 1.3x106
tetrahedral elements. The time step was halved for this computation. The
computed thrust forces are shown in Fig.
12. It can be seen that the agreement between the two analyses is
very good; even the coarse grid produces adequate resolution.

View larger version (26K):
[in this window]
[in a new window]
|
Fig. 12. Effect of grid refinement on the computed thrust. The finer grid had twice
the resolution and the time step was halved.
|
|
Viscous effect
To the study the effects of viscosity, a laminar viscous computation was
carried out, for Re=120. Because of the lack of information on the
transition to and the presence of uncertainties in turbulence modeling,
laminar flow was assumed first. Fig.
13 shows the time history of the thrust and drag forces for the
inviscid and the viscous cases. The finer mesh employed in the grid-refinement
study (see above) was used for this computation, and the mesh size near the
leading and trailing edge of the wing was approximately 0.02. It is clear that
viscous effects are minimal, and the thrust and drag forces are dominated by
translational and rotational mechanisms. Hence, inviscid computations were
carried out to study the effect of phasing between the translational and
rotational motions.
Effect of phasing
Fig. 14 compares the forces
produced when the rotational motion precedes stroke reversal (`advanced' case)
with those of Dickinson et al.
(1999
). Again, the agreement
with the experimental results is good. In this case, the peak in the thrust
force is achieved prior to the beginning of the downstroke at t=10.76
s and is approximately 0.56N compared with a value of 0.47N for the
symmetrical case (see Fig. 3A).
This can be explained by the rotational mechanisms discussed above. The
rotational effect diminishes prior to the beginning of the downstroke,
producing a negative thrust of 0.2N. The thrust then increases until
t=11.98s. During this period, the wing moves through the wake created
during the upstroke, as in the symmetrical case, resulting in a high pressure
on the bottom of the wing. The velocity vectors near the leading edge are
shown in Fig. 15. During this
period, both the translational and rotational accelerations are in phase
(Fig. 8). The peak thrust is
approximately 0.48N compared with a value of 0.28N for the symmetrical case.
Thereafter, the combined effect of rotational and translational motions
produces reduced thrust until a second peak arises due to the rotational
motion at t=14.23s, prior to the beginning of the upstroke. The mean
thrust force for one wingbeat cycle is approximately 0.312N, and the mean
thrust coefficient
is 1.291. The
mean drag force is 0.457N and the drag coefficient
is 1.89.

View larger version (26K):
[in this window]
[in a new window]
|
Fig. 14. Thrust (A) and drag (B) forces for the `advanced' case in which the wing
rotation precedes stroke reversal. The results from the present computational
model are compared with the experimental data of Dickinson et al.
(1999 ).
|
|

View larger version (41K):
[in this window]
[in a new window]
|
Fig. 15. Velocity vectors near the leading edge early in the downstroke for the
`advanced' wing rotation in which rotation precedes stroke reversal. Velocity
vectors are shown in the xz plane at y=10 cm. The vectors
are colored according to the magnitude of velocity (in cm s-1) and
are of constant length. L.E., leading edge.
|
|
In the `delayed' case, where wing rotation is delayed with respect to
stroke reversal, the rotational motion does not produce any thrust prior to
the beginning of the downstroke (Fig.
16A). The mean thrust force for one wingbeat cycle is
approximately 0.206N and the mean thrust coefficient
is 0.854. The
mean drag force is 0.457N and the drag coefficient
is 1.496
(Fig. 16B). In the initial
period following stroke reversal, the rotational effect continues to produce a
negative thrust. Fig. 17A,B
shows the velocity vectors near the leading edge. The leading edge vortex from
the upstroke is not present after t=12.05 s. In this case, the high
pressure on the bottom of the wing together with the orientation of the wing
cause a reduction in thrust. Subsequently, the combined translational and
rotational mechanisms result in an increase in thrust. At t=12.8 s,
we observe a plateau region in the thrust
(Fig. 16A). During this
period, the presence of a trailing edge vortex on the upper surface
(Fig. 17C,D) increases the
pressure on the upper surface of the wing, resulting in a temporary loss of
thrust; when this vortex leaves the trailing edge, thrust increases again.

View larger version (25K):
[in this window]
[in a new window]
|
Fig. 16. Thrust (A) and drag (B) forces for the `delayed' case in which wing
rotation is delayed with respect to stroke reversal. The results from the
present computational model are compared with the experimental data of
Dickinson et al. (1999 ).
|
|

View larger version (72K):
[in this window]
[in a new window]
|
Fig. 17. Velocity vectors near the leading edge (A,B) and the trailing edge (C,D)
early in the downstroke for the `delayed' case of wing motion in which
rotation is delayed relative to stroke reversal. Velocity vectors are shown in
the xz plane at y=10 cm. The vectors are colored according
to the magnitude of velocity (cm s-1) and are of constant length.
L.E., leading edge; T.E., trailing edge.
|
|
Fig. 18 shows the magnitude
of velocity in the xz plane at y=10 cm at the beginning of the
downstroke for the three cases of wing motion in the wake created by the wing.
Velocities are greatest for the advanced case and smallest for delayed
rotation. The wing moving through the higher-velocity fluid therefore produces
an additional thrust in the advanced rotation case, whereas the wing for the
delayed case intercepts the flow at an angle that produces negative thrust.
Similar velocity fields can be seen in the particle image velocimetry data of
Dickinson et al. (1999
).

View larger version (66K):
[in this window]
[in a new window]
|
Fig. 18. Magnitude of velocity in the wake of the wing at the beginning of the
downstroke for the three cases of wing motion: rotation advanced (A),
symmetrical (B) or delayed (C) relative to stroke reversal. Velocity contours
are shown in the xz plane at y=10 cm. The contours are colored
according to the magnitude of absolute velocity (in cm s-1).
|
|
Mechanical aspects of the flapping wing
The results of the present study were then used to derive the input forces,
moments, power requirement and efficiency for the creation of a robotic fly.
First, the forces on the wing were integrated to a particular spanwise
location from the root. The forces were obtained by marking this location and
then computing the force contribution of all the surface elements on the wing
up to this location. This location was then tracked using the prescribed rigid
body motion. The elements contributing to the force up to this location are
recomputed as the surface mesh is regenerated due to the motion. Three
different spanwise locations were chosen: quarter span, half span and
three-quarter span of the wing. Fig.
19 shows their force contributions compared with the total
contribution of the wing and half the total thrust force generated by the wing
for the symmetrical rotation case. It is clear from this figure that nearly
half the thrust is generated by the outer 25% of the wing.

View larger version (28K):
[in this window]
[in a new window]
|
Fig. 19. Spanwise contribution to thrust for the `symmetrical' case of wing rotation
relative to stroke reversal. The thrust generated at three spanwise locations
(quarter, half and three-quarter span) was calculated and is shown together
with the total thrust produced by the wing and half the total thrust.
|
|
Moments about the wing root in the wing coordinate system
(x',y',z') were then computed from the moments
about the fixed body coordinates and the prescribed kinematics of the wing
(Fig. 20). The moment about
the wing rotation axis z' (Mz') is
nearly zero throughout the cycle, implying that there is no torsional load on
the system. The variation of moment about the x' (=x)
axis (Mx') is anti-symmetrical because only one wing
is considered for the moment computation.

View larger version (25K):
[in this window]
[in a new window]
|
Fig. 20. Moments about the wing coordinate system
(x',y',z') for the `symmetrical' case of wing
rotation. Mx', My',
Mz', moment about wing rotation axis
x', y' and z', respectively.
|
|
The power input to the wing Pin is computed by:
 | (12) |
where F is the force vector and Wwing is the velocity of the
surface of the wing. Fig. 21
shows the instantaneous power requirement for one wing for the three phases of
wing rotation. The mean power required is 0.024 W for the symmetrical case,
0.039 W for the advanced case and 0.024 W for the delayed case. The mean
thrust for the symmetrical case is 0.318 N; values for the advanced and
delayed cases are 0.312 N and 0.206 N, respectively.

View larger version (32K):
[in this window]
[in a new window]
|
Fig. 21. Power requirement for one wing for three different phases between wing
rotation and stroke reversal.
|
|
 |
Acknowledgments
|
---|
This work was supported by the Office of Naval Research through the
Tactical Electronic Warfare Division Micro Air Vehicles Program of the Naval
Research Laboratory. The authors would like to thank Professor Michael
Dickinson and Mr Sanjay Sane for providing the experimental kinematics and
data for comparison and for very useful discussions and Professor Rainald
Löhner of the George Mason University for his support throughout the
course of this work. The computations carried out for this work were supported
in part by a grant of HPC time from the DoD HPC centers, ARL MSRC SGI-O2K and
NRL SGI-O2K.
 |
References
|
---|
Anderson, J. M. (1996). Vorticity
control for efficient propulsion. PhD dissertation, Massachusetts
Institute of Technology.
Dickinson, M. H., Lehmann, F.-O. and Sane, S.
(1999). Wing rotation and the aerodynamic basis of insect flight.
Science 284,1954
-1960.[Abstract/Free Full Text]
Ellington, C. P. (1984). The aerodynamics of
hovering insect flight. IV. Aerodynamic mechanisms. Phil. Trans. R.
Soc. Lond. B 305,79
-113.
Freymuth, P. (1999). Thrust generation by an
airfoil in hover modes. Exp. Fluids
9, 17-24.
Jones, K. D. and Platzer, M. F. (1997).
Numerical computation of flappingwing propulsion and power extraction.
AIAA Paper No. 97-0826.
Liu, H. and Kawachi, K. (1998). A numerical
study of insect flight. J. Comp. Phys.
146,124
-156.
Martin, D. and Löhner, R. (1992). An
implicit linelet-based solver for incompressible flows. AIAA
Paper No. 92-0668.
Meirovitch, L. (1970). Methods of
Analytical Dynamics. New York: McGraw-Hill.
Ramamurti, R. and Löhner, R. (1992).
Evaluation of an incompressible flow solver based on simple elements. In
Advances in Finite Element Analysis in Fluid Dynamics,
FED vol. 137 (ed. M. N. Dhaubhadel, M. S. Engleman and
J. S. Reddy), pp. 33-42. New York: ASME
Publication.
Ramamurti, R., Löhner, R. and Sandberg, W. C.
(1994). Evaluation of scalable three-dimensional incompressible
finite element solver. AIAA Paper No.
94-0756.
Ramamurti, R., Löhner, R. and Sandberg, W. C.
(1996). Computation of unsteady flow past a tuna with caudal fin
oscillation. In Advances in Fluid Mechanics, vol.9
(ed. M. Rahman and C. A. Brebbia), pp.169
-178. Southampton: Computational Mechanics
Publications.
Ramamurti, R. and Sandberg, W. C. (2001).
Simulation of flow about flapping airfoils using a finite element
incompressible flow solver. AIAA J.
39,253
-260.
Ramamurti, R., Sandberg, W. C. and Löhner, R.
(1995). Simulation of a torpedo launch using a three-dimensional
incompressible finite element flow solver and adaptive remeshing.
AIAA Paper No. 95-0086.
Walker, J. A. and Westneat, M. W. (1997).
Labriform propulsion in fishes: kinematics of flapping aquatic flight in bird
wrasse Gomphosus varius (Labridae). J. Exp.
Biol. 200,1549
-1569.[Abstract/Free Full Text]