The Journal of Experimental Biology 205, 1507–1518 (2002) 1507

Printed in Great Britain © The Company of Biologists Limited 2002

JEB3871

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: [email protected]

Accepted 4 March 2002

Summary
A finite element flow solver was employed to compute the peak in the thrust forces is higher than when the wing
unsteady flow past a three-dimensional Drosophila wing rotation is in phase with the stroke reversal and that the
undergoing flapping motion. The computed thrust and peak thrust is reduced further when the wing rotation is
drag forces agreed well with results from a previous delayed. As suggested by previous authors, we observe
experimental study. A grid-refinement study was that the rotational mechanism is important and that the
performed to validate the computational results, and a combined translational and rotational mechanisms are
grid-independent solution was achieved. The effect of necessary to describe accurately the force time histories
phasing between the translational and rotational motions and unsteady aerodynamics of flapping wings.
was studied by varying the rotational motion prior to the
stroke reversal. It was observed that, when the wing Key words: unsteady mechanism, incompressible fluid, flapping foil,
rotation is advanced with respect to the stroke reversal, unstructured mesh, insect, flight, Drosophila.

Introduction
Flapping foil propulsion has received considerable attention actual three-dimensional insect wings. Ramamurti et al. (1996)
in the past few years as an alternative to the propeller. This computed three-dimensional unsteady flow past moving and
mode of propulsion, which involves no body undulation, has deforming geometries in a simulation of a swimming tuna with
many applications, such as propulsion of submersibles, caudal fin oscillation.
maneuvering and flow control, that are of interest to the The three-dimensional wing strokes of insects can be
hydrodynamic community and unconventional aerodynamics divided into two translational phases and two rotational phases.
of micro aerial vehicles (MAVs) and the study of aircraft flutter During the translational phases, the upstroke and the
that are of interest to the aerodynamic community. downstroke, the wings move through the air with high angles
Flapping foil propulsion is also important in the area of bio- of attack, and during the rotational phases, pronation and
fluid dynamics for the study of propulsion in insects, birds and supination, the wings rotate rapidly and reverse direction.
certain aquatic animals. Flying animals generate lift and thrust Dickinson et al. (1999) studied the effects of translational and
as a consequence of the interaction between the flapping rotational mechanisms of the wing in Drosophila
motions of the wings and the surrounding air. These animals melanogaster. They directly measured the forces produced by
also perform maneuvers involving rapid plunging and pitching flapping wings and explained the aerodynamics of insect flight
motions. Conventional steady-state theories do not predict by interactions between three unsteady flow mechanisms. The
sufficient forces to meet those required for flight (Ellington, ‘delayed’ stall mechanism is a translational mechanism which,
1984). Therefore, we need to understand the unsteady in two dimensions, produces high lift in the initial phases of
aerodynamics of flapping wings undergoing highly three- translation until eventual flow separation; in three dimensions,
dimensional motions with widely varying geometries. the spanwise flow effectively prevents stall. Rotational
Experimental work on two-dimensional flapping foils has circulation and wake capture are rotational mechanisms that
been carried out by Anderson (1996) and Freymuth (1999). depend mainly on the pronation and supination of the wing
Computational studies have been performed by Jones and during stroke reversal. Walker and Westneat (1997) have
Platzer (1997) and Ramamurti and Sandberg (2001). While studied experimentally the kinematics of fin motion in a fish,
two-dimensional wing section investigations can yield useful the bird wrasse Gomphosus varius. Liu and Kawachi (1998)
insights on the coupled pitching and heaving dynamics, numerically investigated the flow over a hovering hawkmoth,
nothing can be learned concerning the influence of spanwise Manduca sexta. They reported the presence of a spiral leading-
flow. It is therefore essential to carry out computations for edge vortex to which they attributed the lift enhancement
1508 R. Ramamurti and W. C. Sandberg
mechanism. They validated their results by comparing the A Thrust (lift)
computational streak lines found in a two-dimensional
hovering airfoil with the experimental smoke visualization, but
Upstroke
they did not directly compare instantaneous forces.
Here, we extend the two-dimensional pitching and heaving Stroke plane
y
airfoil computations to three dimensions. This study will Downstroke
address the role of the rotational motion in detail. Also, the role z Drag
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 x
computations by comparing the forces with the experimental
results of Dickinson et al. (1999) and performing grid-
refinement studies, to verify the hypothesis of Dickinson et al. B
(1999) that rotational mechanisms of the wing form the basis y Beginning of downstroke
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. x, x'
To this end, computations are performed for various phase y'
angles between the rotation and translation motions, and the
z Mid-stroke
time history of the unsteady forces is compared with the
φ
experimental results. The flow solver we employ is a finite-
element-based incompressible flow solver based on simple,
R
low-order elements. The simple elements enable the flow solver
z'
to be as fast as possible, reducing the overheads in building
End of downstroke
element matrices, residual vectors, etc. The governing
equations are written in Arbitrary Lagrangian Eulerian form,
Fig. 1. (A) Schematic diagram of a hovering Drosophila showing the
which enables flow with moving bodies to be simulated. The
orientation of the x,y,z coordinate system. (B) Schematic diagram of
details of the flow solver, the rigid body motion and adaptive the flapping Drosophila wing. The position of the wing is shown at
remeshing are given by Ramamurti et al. (1995) and are three different times during the flapping cycle. The coordinate
summarized below. 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.
Materials and methods
The incompressible flow solver discretized in space using a Galerkin procedure with linear
The governing equations employed are the incompressible tetrahedral elements. To be as fast as possible, the overheads
Navier–Stokes equations in Arbitrary Lagrangian Eulerian in building element matrices, residual vectors, etc., should be
(ALE) formulation. They are written as: kept to a minimum. This requirement is met by employing
simple, low-order elements that have all the variables (v,p) at
dv
+ va · ∇v + ∇p = ∇ · σ , (1) the same location. The resulting matrix systems are solved
dt iteratively using a preconditioned conjugate gradient
algorithm (PCG), as described by Martin and Löhner (1992).
dv ∂v
= + w · ∇v , (2) The flow solver has been successfully evaluated for both two-
dt ∂t dimensional and three-dimensional laminar and turbulent flow
problems by Ramamurti and Löhner (1992) and Ramamurti et
∇·v=0. (3)
al. (1994).
Here, p denotes pressure, and va=v–w, the advective velocity To carry out computations of the flow about oscillating and
vector, where v is flow velocity and w is mesh velocity and deforming geometries, one needs to describe grid motion on
the material derivative is with respect to the mesh velocity w. a moving surface, i.e. to couple the moving surface grid to the
Both the pressure p and the stress tensor σ have been volume grid. The volume grid in the proximity of the moving
normalized by the (constant) density ρ and are discretized in surface is then remeshed to eliminate badly distorted
time t using an implicit time-stepping procedure. It is elements. The velocity of the mesh is obtained in a manner so
important for the flow solver to be able to capture the as to reduce this distortion. A detailed description of the
unsteadiness of a flow field. The present flow solver is built various mesh movement algorithms is given in Ramamurti et
as time-accurate from the onset, allowing local time stepping al. (1994). In that study, smoothing of the coordinates was
as an option. The resulting expressions are subsequently employed for the mesh movement with a specified number of
 Three-dimensional computational study of insect flight 1509
Pitch angle 0.8
Advanced Present study A
Roll angle φ Symmetrical
Dickinson et al. (1999)
Delayed
100 60 0.6 0 4 8
A 3 7
40
0.4

Thrust (N)
50

Pitch angle (degrees)

Roll angle (degrees)

20
0.2
0 0

2 6
–20 0
1 5 9
–50
Downstroke Upstroke
–40
–0.2
10 12 14 16 18 20
Downstroke Upstroke
–100 –60
10 12 14 16 18 20 1.5
Present study B
Angular velocity Dickinson et al. (1999)
Advanced
Translational velocity 1.0
Symmetrical
Delayed
25.0 100
B 0.5
Drag (N)
Translational velocity (cm s–1)

12.5 50
Angular velocity (° s–1)

0 0 –0.5

Downstroke Upstroke
–12.5 –50 –1.0
10 12 14 16 18 20
Time (s)
Downstroke Upstroke Fig. 3. Time history of thrust (A) and drag (B) forces during one
–25.0 –100
10 12 14 16 18 20 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
Time (s)
to times discussed in the text.
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 To reduce the distortion of the mesh, the coordinates at the
velocity of the wing tip and rotational (angular) velocity of the wing new time xn+1 were obtained as a weighted average of the
for three different phases. original grid point location at time t=0(x0) and the location of
the point as if it moved rigidly with the body (xn+1
rigid):

layers of elements that move rigidly with the wing. In two- xn+1 = x0f(r) + xn+1
rigid[1 − f(r)] , (4)
dimensional studies (Ramamurti and Sandberg, 2001), the where the weighting function (f ) is a simple linear function
grid showed that the elements at the edge of the rigid layers based on the distance from the center of rotation r, and is given
were very distorted after one cycle of oscillation. This is due by:
to a residual mesh velocity that is present as a result of the
non-convergence of the mesh velocity field. This will appear f(r) = 0 for r < rmin , (5)
whether a spring-analogy or a Laplacian-based smoothing is
used. f(r) = 1 for r > rmax (6)
1510 R. Ramamurti and W. C. Sandberg
y
A

Leading Trailing
edge edge

y=10 cm

x
Wing chord

z
Downstroke
Fig. 5. Instantaneous particle traces at the beginning of the
downstroke. Particles were released from a rake of rectangular grid
y
of points in a plane 0.8 mm away from the bottom surface of the
B x=0 cm
wing. Using the instantaneous velocity field, the positions of these
particles were obtained by integrating the velocity at these rake
Tip points until the length of the traces exceeded a specified length or the
particles ended on a solid boundary or exited the computational
Leading Trailing domain. These particle traces are colored according to the magnitude
edge edge of the velocity (in cm s–1) at that location. A leading edge vortex is
y=10 cm seen rotating in the counterclockwise direction, and a stagnation line
is shown near the z′ axis of rotation (dark blue traces).
x
Wing chord pulled away from the symmetry plane. To avoid this problem,
the points on the symmetry plane are allowed to glide along
z this plane. A similar technique has been employed for the
Downstroke simulation of torpedo launch from a submarine by Ramamurti
et al. (1995) where the gap between the launch tube and the
torpedo was small.
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 Results and discussion
plane for which the velocity vectors are shown in Fig. 6 is indicated.
The configuration of the hovering Drosophila
(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.
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
and downstroke), the wing moves from close to the y axis through
(r − rmin) an angle φ, the wingbeat amplitude. The flapping wing
f(r) = for rmax > r > rmin . (7)
(rmax − rmin) configuration used in the flow simulations is shown in
Fig. 1B. This is based on the experimental arrangement of
The mesh velocity is then obtained using: Dickinson et al. (1999). Fig. 1B shows the position of the
1 wing at three different times during the flapping cycle. The
w= (xn=1 − xn) . (8) coordinate system (x′,y′,z′) is fixed to the wing, and the wing
∆t rotates about the z′ axis throughout the cycle. The planform
The values of rmin and rmax used in this study are 2.0 and of the wing is derived from the Drosophila wing and is 25 cm
10.0, respectively. To reduce the computational effort, the long and 3.2 mm thick. The experimental apparatus consisted
experimental arrangement of Dickinson et al. (1999) was of two wings immersed in a tank of mineral oil. The viscosity
approximated by introducing a symmetry plane. Because of the of the oil, the length of the wing and the frequency of the
proximity of the wing at the beginning of the downstroke and flapping motion were chosen to match the Reynolds number
the rotation of the wing during the pronation phase, the normal (Re) of a typical Drosophila melanogaster, approximately
component mesh velocity of the points on the symmetry plane 136. The Re for the present calculations is defined on the
can become non-zero. This would result in the points being basis of the mean chord of the wing c– (6.7 cm) and the mean
 Three-dimensional computational study of insect flight 1511

A B C

x x
L.E. x L.E. L.E.

z z
z

t=11.81 s t=11.897 s t=11.983 s

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.

A B C

T.E. T.E.
T.E.

x x x

z
z z
t=12.24 s t=12.5 s t=12.672 s

C E F

Wing tip
Wing tip
Wing tip

y y y

x x x

t=12.24 s t=12.5 s t=12.672 s

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.
1512 R. Ramamurti and W. C. Sandberg
Angular acceleration
A
Advanced
Translational Symmetrical
acceleration Delayed
100 200

150
Translational acceleration (cm s–2)

Angular acceleration (° s–2)

50 100
3
50

0 0

–50

–50 –100

1 –150 B
Downstroke Upstroke
–100 –200
10 12 14 16 18 20
Time (s)
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.

wing-tip velocity Ut (ignoring the forward velocity), as

follows:
c–Ut
Re = , (9)
ν
Fig. 9. Particle streaklines produced during the downstroke at
where c–=2R/AR, Ut=2φnR, R is the wing length (25 cm), AR is t=12.93 s. The particles were release from a rake of rectangular grid
the aspect ratio of the wing, n is the frequency of flapping on a plane 3.0 cm above (A) and 3.5 cm below (B) the wing parallel
motion, φ is the wingbeat amplitude (peak to peak, in rad) and to the wing and near the leading edge. The absolute velocities of the
ν is the kinematic viscosity (115 cSt=115×10–6 m2 s–1). particles (in cm s–1) are given in the color scale.
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 and rotation. First, an inviscid solution was obtained
Fig. 2A. The motion of the wing is prescribed using these two using a grid consisting of 178 219 points and 965 877
angles. Fig. 2B shows the translational velocity of the wing tetrahedral elements. An initial steady-state solution was
tip and the rotational (angular) velocity of the wing. Three obtained in 1500 time steps. The unsteady solution using the
different phases between the translational and rotational prescribed kinematics (Fig. 2) is then obtained. The surface
motions were used. In the ‘advanced’ case, wing rotation pressure on the wing is integrated to obtain the forces on the
precedes stroke reversal by 8 % of the wingbeat cycle; wing along the three axes (Fx, Fy, Fz). The thrust T and the
‘symmetrical’ wing motion is where the wing rotation drag D forces are then computed as T=–Fx and
occurs symmetrically with respect to stroke reversal; in D=√(Fy2+Fz2), respectively. These forces are compared with
‘delayed’ wing motion, rotation is delayed by 8 % with those obtained from the experiments of Dickinson et al.
respect to stroke reversal. Wingbeat amplitude is 160 °, (1999).
flapping frequency is 0.145 Hz and the angle of attack at The unsteady computation was carried out for five cycles of
midstroke is approximately 40 ° during both upstroke and oscillation. Fig. 3 shows the thrust and drag forces during one
downstroke. cycle of the wingbeat and compares the values with those of
Dickinson et al. (1999). The present computations capture the
Symmetrical case peak forces well. The mean thrust force is approximately
–
The flow solver described here is employed to compute 0.318 N, and the mean thrust coefficient Cr is 1.317. The mean
–
the flow past the Drosophila wing undergoing translation drag force is 0.375 N and the drag coefficient Cp is 1.55. The
 Three-dimensional computational study of insect flight 1513

A B

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

force coefficients were obtained using the following non- If a rotational mechanism alone were present, the thrust should
dimensionalization: continue to decrease until t3; in fact, the thrust force increases
– between t1 and t2. This happens after the wing changes direction
T
CT = , (10) at the start of each half-stroke. Dickinson et al. (1999) attribute
1
2– 2 this increase in thrust to a wake-capture mechanism in which the
ρUt c Rr2(s)
2 wing passes through the shed vorticity of the previous stroke.
– The position of the wing at the beginning of the downstroke
D
CD = , (11) is shown in Fig. 4A. The chord in this case of symmetrical
1 rotation is aligned with the x axis at this instant. We found a
2– 2
ρUt c Rr2(s)
2
– –
where T and D are the mean thrust and drag forces, respectively,
2
and r2 (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 Fig. 11. Particle traces prior to the end of the downstroke, t=14.65 s.
continuously rotates in a counterclockwise direction producing a The particles were release from a rake of rectangular grid on a plane
circulation pointing nearly along the +y direction. Between t0 and 3.2 cm above the wing parallel to the wing and near the leading edge.
t1, the wing is translating in the –z direction, resulting in a force Traces are colored according to the magnitude of absolute velocity
pointing in the –x direction, thus producing a peak in thrust at t0. (in cm s–1).
1514 R. Ramamurti and W. C. Sandberg
0.5 0.8
Present study A
0.4 Dickinson et al. (1999)
0.6

0.3
0.4
Thrust (N)

Thrust (N)
0.2
0.2
0.1
0
0
Finer grid
Coarse grid Downstroke Upstroke
–0.1 –0.2
10 12 14 16 18 20 10 12 14 16 18 20
Time (s)
1.5
Fig. 12. Effect of grid refinement on the computed thrust. The finer Present study
grid had twice the resolution and the time step was halved. Dickinson et al. (1999)
B
1.0

0.8
Inviscid A 0.5

Drag (N)
Viscous
0.6
0

0.4
Thrust (N)

–0.5

0.2 Downstroke Upstroke

–1.0
10 12 14 16 18 20
0 Time (s)

Downstroke Upstroke Fig. 14. Thrust (A) and drag (B) forces for the ‘advanced’ case in
–0.2 which the wing rotation precedes stroke reversal. The results from
10 12 14 16 18 20 the present computational model are compared with the experimental
data of Dickinson et al. (1999).
1.5
Inviscid B
Viscous
1.0 directions until the length of the traces exceeded a specified
length or the particles ended on a solid boundary or exited the
0.5 computational domain. These instantaneous particle traces are
Drag (N)

colored according to the magnitude of velocity at that location.

Fig. 5 shows the leading edge vortex with vorticity oriented in
0
the +y direction. This leading edge vortex was created at the end
of the preceding upstroke. This vortex is located below the wing
–0.5 and is rotating in the counterclockwise direction, which can be
Downstroke Upstroke seen from the velocity vectors shown in Fig. 6A. A possible
–1.0 explanation for the increase in thrust between t1 and t2 is that the
10 12 14 16 18 20 wing moving through this wake benefits from the shed vorticity.
Time (s) As the wing moves through this vortex during the downstroke, it
Fig. 13. Effects of viscosity on the computed thrust (A) and drag (B) produces a stagnation region at the bottom of the wing near the
forces, Re=120. z′ axis (Fig. 5), resulting in an increase in thrust.
At the beginning of the downstroke (t1), the flow separates at
the leading edge and reattaches on the bottom of the wing as
separation bubble attached to the leading edge during the shown in Fig. 6A. As the wing continues to move down, between
interval t0–t1. Particles were released from a rake of rectangular t1 and t2, the separation point of this bubble moves back along
grid of points in a plane 0.8 mm away from the bottom surface the wing chord (Fig. 6). During the interval t2–t3, we observe a
of the wing. Using the instantaneous velocity field, the positions trailing edge separation bubble forming (see Fig. 7A–C). A
of these particles were obtained by integrating the velocity at similar separation region forms at the wing tip, as can be seen in
these rake points in both the positive and negative velocity Fig. 7D–F. Fig. 4B shows the position of the wing at t=12.5 s.
 Three-dimensional computational study of insect flight 1515

A B C

x
x x

L.E. L.E.
L.E.
z z
z

t=11.638 s t=11.81 s t=11.983 s

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.

Thrust production reaches a local minimum around this instant 0.8

(Fig. 3). In Fig. 7A,D, a large recirculation region can be seen in
0.6
Present study A
the wake of the wing. This separated flow from the wing tip and Dickinson et al. (1999)
trailing edge will result in a higher pressure on the upper surface 0.4
of the wing and, hence, a reduction in the thrust. During the 0.2
interval t1–t3, the magnitude of the translational acceleration of
Thrust (N)

the wing decreases while that of the angular acceleration 0

increases (symmetrical phase, Fig. 8). –0.2
Between t2 and t3, the magnitude of the translational
acceleration is large enough to overcome the rotational effect –0.4
and, when the angular acceleration becomes large enough, the –0.6
rotational mechanism dominates, resulting in the observed Downstroke Upstroke
–0.8
reduction in thrust until t3. Between t3 and t4, the translational 10 12 14 16 18 20
effect should result in a constant thrust because the
translational acceleration is almost constant during this period. 1.5
The rotational effect produces an increase in thrust between t3 B
and t4, with a plateau in the middle, which occurs when the 1.0
trailing edge vortex is shed. Similar trends are observed during
the supination phase prior to the beginning of the upstroke
0.5
(t4–t5) and at the beginning of the upstroke (t5–t7).
Drag (N)

Fig. 9 shows the instantaneous traces or streaklines of

particles released 3.0 cm above or 3.5 cm below the wing in a 0
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 –0.5 Present study
Dickinson et al. (1999)
visible above the wing surface. A leading edge vortex spinning
Downstroke Upstroke
in the counterclockwise direction is found below the wing –1.0
surface (Fig. 9B). Particle traces near the mid stroke are shown 10 12 14 16 18 20
in Fig. 10A. At this instant, the wing rotation axis z′ is aligned Time (s)
with the body coordinate z (see Fig. 1B). Here, we can see the Fig. 16. Thrust (A) and drag (B) forces for the ‘delayed’ case in
beginnings of the leading edge vortex on the upper surface of which wing rotation is delayed with respect to stroke reversal. The
the wing that is also shown by the velocity vectors in the xy plane results from the present computational model are compared with the
at z=20 cm (Fig. 10B). Fig. 10C shows the pressure contours on experimental data of Dickinson et al. (1999).
the upper surface of the wing. The pressure is non-
dimensionalized with respect to the dynamic head, GρUt2, where
ρ is the density of the mineral oil (880 kg m–3) used in the pressure is present from the root of the wing up to approximately
experiments of Dickinson et al. (1999). A region of constant 60 % of the span and extending to the trailing edge. Fig. 11
1516 R. Ramamurti and W. C. Sandberg
shows the particle traces prior to the end of the downstroke. can be explained by the rotational mechanisms discussed above.
Here, the leading edge vortex is clearly seen on the upper surface The rotational effect diminishes prior to the beginning of the
of the wing. downstroke, producing a negative thrust of 0.2 N. The thrust then
increases until t=11.98 s. During this period, the wing moves
Grid-refinement study through the wake created during the upstroke, as in the
To assess the sensitivity of our computational results, we symmetrical case, resulting in a high pressure on the bottom of
carried out a grid-refinement study. The resolution of the grid the wing. The velocity vectors near the leading edge are shown
in the vicinity of the wing was doubled. The computations were in Fig. 15. During this period, both the translational and rotational
carried out using a grid consisting of approximately 238×103 accelerations are in phase (Fig. 8). The peak thrust is
points and 1.3×106 tetrahedral elements. The time step was approximately 0.48 N compared with a value of 0.28 N for the
halved for this computation. The computed thrust forces are symmetrical case. Thereafter, the combined effect of rotational
shown in Fig. 12. It can be seen that the agreement between and translational motions produces reduced thrust until a second
the two analyses is very good; even the coarse grid produces peak arises due to the rotational motion at t=14.23 s, prior to the
adequate resolution. beginning of the upstroke. The mean thrust force for one
wingbeat cycle is approximately 0.312 N, and the mean thrust
Viscous effect –
coefficient CT is 1.291. The mean drag force is 0.457 N and the
–
To the study the effects of viscosity, a laminar viscous drag coefficient CD is 1.89.
computation was carried out, for Re=120. Because of the lack In the ‘delayed’ case, where wing rotation is delayed with
of information on the transition to and the presence of respect to stroke reversal, the rotational motion does not
uncertainties in turbulence modeling, laminar flow was assumed produce any thrust prior to the beginning of the downstroke
first. Fig. 13 shows the time history of the thrust and drag forces (Fig. 16A). The mean thrust force for one wingbeat cycle is
–
for the inviscid and the viscous cases. The finer mesh employed approximately 0.206 N and the mean thrust coefficient Cr is
in the grid-refinement study (see above) was used for this 0.854. The mean drag force is 0.457 N and the drag coefficient
–
computation, and the mesh size near the leading and trailing CD is 1.496 (Fig. 16B). In the initial period following stroke
edge of the wing was approximately
0.02. It is clear that viscous effects are A B
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.
L.E.
Effect of phasing
Fig. 14 compares the forces produced L.E.
when the rotational motion precedes x
x
stroke reversal (‘advanced’ case) with
those of Dickinson et al. (1999). Again,
z
the agreement with the experimental z
results is good. In this case, the peak in
the thrust force is achieved prior to the t=11.81 s t=12.07 s
beginning of the downstroke at t=10.76 s C D
and is approximately 0.56 N compared
with a value of 0.47 N for the
symmetrical case (see Fig. 3A). This

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 T.E.
T.E.
delayed relative to stroke reversal. Velocity x
vectors are shown in the xz plane at x
y=10 cm. The vectors are colored
according to the magnitude of velocity
z
(cm s–1) and are of constant length. L.E., z
leading edge; T.E., trailing edge. t=12.59 s t=12.76 s
 Three-dimensional computational study of insect flight 1517

A B

Advanced Symmetrical

Fig. 18. Magnitude of velocity in the wake of the wing at the

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

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).
Delayed 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
0.8
this vortex leaves the trailing edge, thrust increases again.
Thrust (full span) Fig. 18 shows the magnitude of velocity in the xz plane at
Half thrust
3/4 span y=10 cm at the beginning of the downstroke for the three cases
0.6 1/2 span
1/4 span
0.10
0.4 Mx′
Thrust (N)

My′
0.05 Mz′
0.2
M (N m)

0
0

-0.05
–0.2
10 12 14 16 18 20 Downstroke Upstroke
Time (s) –0.10
10 12 14 16 18 20
Fig. 19. Spanwise contribution to thrust for the ‘symmetrical’ case of Time (s)
wing rotation relative to stroke reversal. The thrust generated at three
spanwise locations (quarter, half and three-quarter span) was Fig. 20. Moments about the wing coordinate system (x′,y′,z′) for the
calculated and is shown together with the total thrust produced by the ‘symmetrical’ case of wing rotation. Mx′, My′, Mz′, moment about
wing and half the total thrust. wing rotation axis x′, y′ and z′, respectively.
1518 R. Ramamurti and W. C. Sandberg
0.10 The power input to the wing Pin is computed by:
Symmetrical
Advanced dφ dα
0.08 Delayed Pin = F · wwing = Mx′ + Mz′ , (12)
dt ∂t
where F is the force vector and wwing is the velocity of the
0.06
surface of the wing. Fig. 21 shows the instantaneous power
Input power (W)

requirement for one wing for the three phases of wing rotation.
0.04 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
0.02 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.
0
Downstroke Upstroke This work was supported by the Office of Naval Research
–0.02 through the Tactical Electronic Warfare Division Micro Air
10 12 14 16 18 20 Vehicles Program of the Naval Research Laboratory. The
Time (s) authors would like to thank Professor Michael Dickinson and
Mr Sanjay Sane for providing the experimental kinematics
Fig. 21. Power requirement for one wing for three different phases
between wing rotation and stroke reversal.
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
of wing motion in the wake created by the wing. Velocities are computations carried out for this work were supported in part
greatest for the advanced case and smallest for delayed by a grant of HPC time from the DoD HPC centers, ARL
rotation. The wing moving through the higher-velocity fluid MSRC SGI-O2K and NRL SGI-O2K.
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 References
et al. (1999). 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
Mechanical aspects of the flapping wing the aerodynamic basis of insect flight. Science 284, 1954–1960.
The results of the present study were then used to derive the Ellington, C. P. (1984). The aerodynamics of hovering insect flight. IV.
Aerodynamic mechanisms. Phil. Trans. R. Soc. Lond. B 305, 79–113.
input forces, moments, power requirement and efficiency for Freymuth, P. (1999). Thrust generation by an airfoil in hover modes. Exp.
the creation of a robotic fly. First, the forces on the wing were Fluids 9, 17–24.
integrated to a particular spanwise location from the root. The Jones, K. D. and Platzer, M. F. (1997). Numerical computation of flapping-
wing propulsion and power extraction. AIAA Paper No. 97-0826.
forces were obtained by marking this location and then Liu, H. and Kawachi, K. (1998). A numerical study of insect flight. J. Comp.
computing the force contribution of all the surface elements on Phys. 146, 124–156.
the wing up to this location. This location was then tracked Martin, D. and Löhner, R. (1992). An implicit linelet-based solver for
incompressible flows. AIAA Paper No. 92-0668.
using the prescribed rigid body motion. The elements Meirovitch, L. (1970). Methods of Analytical Dynamics. New York: McGraw-
contributing to the force up to this location are recomputed as Hill.
the surface mesh is regenerated due to the motion. Three Ramamurti, R. and Löhner, R. (1992). Evaluation of an incompressible flow
solver based on simple elements. In Advances in Finite Element Analysis in
different spanwise locations were chosen: quarter span, half Fluid Dynamics, FED vol. 137 (ed. M. N. Dhaubhadel, M. S. Engleman and
span and three-quarter span of the wing. Fig. 19 shows their J. S. Reddy), pp. 33–42. New York: ASME Publication.
force contributions compared with the total contribution of the Ramamurti, R., Löhner, R. and Sandberg, W. C. (1994). Evaluation of
scalable three-dimensional incompressible finite element solver. AIAA
wing and half the total thrust force generated by the wing for Paper No. 94-0756.
the symmetrical rotation case. It is clear from this figure that Ramamurti, R., Löhner, R. and Sandberg, W. C. (1996). Computation of
nearly half the thrust is generated by the outer 25 % of the wing. 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.
Moments about the wing root in the wing coordinate system Southampton: Computational Mechanics Publications.
(x′,y′,z′) were then computed from the moments about the fixed Ramamurti, R. and Sandberg, W. C. (2001). Simulation of flow about
body coordinates and the prescribed kinematics of the wing flapping airfoils using a finite element incompressible flow solver. AIAA J.
39, 253–260.
(Fig. 20). The moment about the wing rotation axis z′ (Mz′) is Ramamurti, R., Sandberg, W. C. and Löhner, R. (1995). Simulation of a
nearly zero throughout the cycle, implying that there is no torpedo launch using a three-dimensional incompressible finite element flow
torsional load on the system. The variation of moment about solver and adaptive remeshing. AIAA Paper No. 95-0086.
Walker, J. A. and Westneat, M. W. (1997). Labriform propulsion in fishes:
the x′ (=x) axis (Mx′) is anti-symmetrical because only one kinematics of flapping aquatic flight in bird wrasse Gomphosus varius
wing is considered for the moment computation. (Labridae). J. Exp. Biol. 200, 1549–1569.