## SUMMARY

We compare computational, experimental and quasi-steady forces in a generic hovering wing undergoing sinusoidal motion along a horizontal stroke plane. In particular, we investigate unsteady effects and compare two-dimensional (2D) computations and three-dimensional (3D) experiments in several qualitatively different kinematic patterns. In all cases, the computed drag compares well with the experiments. The computed lift agrees in the cases in which the sinusoidal changes in angle of attack are symmetrical or advanced with respect to stroke positions, but lags behind the measured 3D lift in the delayed case.

In the range of amplitudes studied here, 3–5 chords, the force coefficients have a weak dependence on stroke amplitude. As expected, the forces are sensitive to the phase between stroke angle and angle of attack, a result that can be explained by the orientation of the wing at reversal. This dependence on amplitude and phase suggests a simple maneuver strategy that could be used by a flapping wing device.

In all cases the unsteady forces quickly reach an almost periodic state with continuous flapping. The fluid forces are dominated by the pressure contribution. The force component directly proportional to the linear acceleration is smaller by a factor proportional to the ratio of wing thickness and stroke amplitude; its net contribution is zero in hovering. The ratio of wing inertia and fluid force is proportional to the product of the ratio of wing and fluid density and the ratio of wing thickness and stroke amplitude; it is negligible in the robotic wing experiment, but need not be in insect flight.

To identify unsteady effects associated with wing acceleration, and coupling between rotation and translation, as well as wake capture, we examine the difference between the unsteady forces and the estimates based on translational velocities, and compare them against the estimate of the coupling between rotation and translation, which have simple analytic forms for sinusoidal motions. The agreement and disagreement between the computed forces and experiments offer further insight into when the 3D effects are important.

A main difference between a 3D revolving wing and a 2D translating wing is
the absence of vortex shedding by a revolving wing over a distance much longer
than the typical stroke length of insects. No doubt such a difference in
shedding dynamics is responsible in part for the differences in *steady
state* force coefficients measured in 2D and 3D. On the other hand, it is
unclear whether such differences would have a significant effect on transient
force coefficients before the onset of shedding. While the 2D steady state
force coefficients underpredict 3D forces, the transient 2D forces measured
prior to shedding are much closer to the 3D forces. In the cases studied here,
the chord is moving between 3 to 5 chords, typical of hovering insect stroke
length, and the flow does not appear to separate during each stroke in the
cases of advanced and symmetrical rotation. In these cases, the wing reverses
before the leading edge vortex would have time to separate even in 2D. This
suggests that the time scale for flow separation in these strokes is dictated
by the flapping frequency, which is dimensionally independent. In such cases,
the 2D unsteady forces turn out to be good approximations of 3D
experiments.

- insect flight
- computational fluid dynamics
- biofluid dynamics
- vorticity field
- two-dimensional force
- three-dimensional force

## Introduction

The demonstrated importance of unsteady effects in insect flight has prompted recent development of better experimental and computational tools to investigate unsteady forces and vortical flows around a flapping wing. In particular, dynamically scaled robots of both the hawkmoth (Ellington et al., 1996) and fruitfly (Dickinson et al., 1999) have been developed to perform controlled measurements of flows and forces. In parallel, researchers have developed direct numerical simulations to probe the detailed vortex dynamics and unsteady forces in flapping flight (Liu et al., 1998; Wang, 2000a,b; Sun and Tang, 2002; Ramamurti and Sandberg, 2002).

Given the complexity of modeling fluid flows in three dimensions (Liu et al., 1998; Sun and Tang, 2002; Ramamurti and Sandberg, 2002), it would be desirable to determine if simpler models provide results that are consistent with those generated in experiments. Here, we compare two-dimensional (2D) computations of hovering flight against robotic wing experiments. 2D computations are appealing partly because of their relative simplicity and efficiency. Obviously, 2D computations cannot predict three-dimensional (3D) effects; on the other hand, it is almost impossible to attach the significance of 3D effects without knowing what happens in 2D flow. Therefore, in addition to being relevant to cases where the flow is approximately 2D, as with large wing aspect ratio, when compared with 3D experiments or computations, 2D computations can offer useful insight into the relative significance of 3D effects, as we will discuss at the end of the paper.

Comparing computations and experiments is delicate, partly because it is almost impossible to match any two setups exactly, and partly because it is tempting to present results that compare well, thereby biasing the interpretations. Therefore, it is essential to test the methods in qualitatively different flows generated by different wing kinematics.

In addition to comparing the experimental and computational forces, we also evaluate the relative importance of unsteady effects. These include wing acceleration, both in translation and rotation, and interactions between the wing the existing flow. Most recent work using a robotic fruit fly focused on kinematics based on tethered flight measurements. These kinematics have relatively constant translational velocity in the mid-stroke and large accelerations and sharp rotations at the end of strokes. In these strokes, pronounced peaks appear near the end of each stroke. These peaks were attributed to either wing rotation and wake capture (Dickinson et al., 1999; Birch and Dickinson, 2003), or rotational and translational acceleration (Sun and Tang, 2002; Sane and Dickinson, 2002). The sinusoidally varying strokes studied here offer a set of kinematics where the relative contribution of some of the dynamic effects can be theoretically estimated. For example, we estimate the relative contribution that wing rotation and acceleration make to the quasi-steady forces. We also estimate the wing inertia relative to the fluid force, as well as the non-inertial forces due to wing acceleration relative to the pressure forces associated with vorticity flux. Given that the free flight kinematics of fruit flies appear to be more sinusoidal than those derived from tethered flight (Fry et al., 2003), our results, though using an idealized kinematics, may nevertheless relate to the forces generated by the real flies.

## Materials and methods

### Experimental method

We obtained both force and flow data using a dynamically scaled robotic
fly, as described in prior studies
(Dickinson and Götz, 1993;
Dickinson et al., 1999;
Birch and Dickinson, 2003). At
the base of one arm was attached a 2D force sensor that measured forces
parallel and perpendicular to the wing surface. The force sensor at the wing
base was designed to be insensitive to the force moments; thus the force
distribution on the wing did not influence overall force magnitude. Lift and
drag forces were then calculated from the perpendicular shear forces measured
by the sensor. The wing consisted of a 2.25 mm thick piece of Plexiglas, cut
to the planform of a *Drosophila* wing. When attached to the force
sensor, total wing length was 0.25 m. The wing and arm apparatus were lowered
into a 1 m×1 m×2 m Plexiglas tank filled with 1.8 m^{3} of
mineral oil (density 0.88×10^{3} kg m^{–3};
kinematic viscosity=115 cSt). By changing flapping frequency we could operate
the robot at Reynolds numbers *Re* between 50 and 200. The design of
the robotic arms permitted the wing to move with three degrees of freedom
(measured as stroke amplitude, stroke deviation from horizontal, and wing
rotation) and a custom program written in Matlab (Mathworks v.5.3) allowed us
to program arbitrary kinematic patterns. To simplify the comparisons with the
2D computational fluid dynamics model, we created a simple back-and-forth wing
beat pattern with no stroke plane deviation.

We used digital particle image velocimetry (DPIV) to measure the flow
structure in a 841 cm^{2} area centered on the wing. The oil was
seeded with air forced through a ceramic water filter stone, creating a dense
bubble field. After the larger bubbles rose to the surface, the remaining
bubbles, although slightly positively buoyant, were effectively stationary for
the duration of each exposure pair. A commercial software package controlling
a dual Nd-YAG laser system (Insight v.3.2, TSI Inc., St Paul, MN, USA) created
two identically positioned light sheets approximately 2.5 mm thick separated
by 2000 μs. These light sheets were parallel to the wing chord and
positioned at 0.65*R*, where *R* is the wing span, and timed to
fire when the wing chord was directly in front of the high-speed video camera
placed perpendicular to the laser sheet. We chose 0.65*R* as our point
of measurement because in a prior DPIV study in which the wake was viewed from
the rear, 0.65*R* was the position in which the circulation was the
greatest (Birch and Dickinson,
2003). We captured one image per stroke from a 29 cm×29 cm
area centered on the wing during each of the four strokes. After saving the
captured images, the trigger for the laser was advanced and the starting
position of the wing was adjusted to line up with the camera at the
appropriate time before starting the next trial
(Birch and Dickinson, 2003). We
repeated this procedure until we had divided each stroke into 10 equally
spaced intervals. In this way, we quantified the fluid flow from the
perspective of the wing through four downstroke/upstroke cycles, although this
paper will only report on the wake dynamics during the fourth stroke. For each
image-pair captured, a 2-frame cross-correlation of pixel intensity peaks with
50% overlap of 64 pixel×64 pixel interrogation areas resulted in 900
velocity vectors/image. Vector validation resulted in the removal of only 2 of
9000 vector values; these were filled by interpolation of a mean value from a
3 × 3 nearest neighbor matrix. A custom program written in Matlab
(v.5.0) calculated vorticity, ω=Δ×*u*, from velocity
fields smoothed using using a least-squares finite difference scheme.

### Computational method

The computation models a thin wing element of elliptic cross section under the same kinematics as performed in the experiments. The computation of flows around this hovering wing employs a fourth-order finite difference scheme of Navier–Stokes equation in vorticity-stream function formulation (E and Liu, 1996). The scheme is implemented in the elliptic coordinates with appropriate boundary conditions to account for the wing motion (Wang, 2000a,b). See also Russell and Wang (2003) for an alternative method employing Cartitian grids appropriate for multiple wings.

The two-dimensional Navier-Stokes equation governing the vorticity in the
elliptic coordinates has the following form:
1
2
where **u** is the velocity field, ω the vorticity field, *v*
is the kinematic viscosity and *S* is the local scaling factor,
*S*(μ,θ)= cosh^{2}μ–cos^{2}θ
resulting from coordinates transformation. The derivatives are with respect to
the elliptic coordinates (μ,θ), whose mesh points are naturally
clustered around the tips and the body of the ellipse to resolve the boundary
layer. It is convenient to express both the vorticity and velocity in terms of
the stream function, Ψ: **u**=–Δ×Ψ andω
=Δ^{2}Ψ. The conformal transformation results in a
constant coefficient Poisson equation for the stream function Ψ, which can
be solved *via* Fast Fourier Transform (FFT) efficiently.

The Navier–Stokes equation is solved in a frame fixed to the wing. In
the 2D vorticity stream function formulation, the non-inertial frame
introduces only one extra term, the rotational acceleration of the wing. Other
non-inertial terms can be expressed as a gradient of a potential function.
Thus they can be absorbed into the pressure term. The curl of the gradient of
pressure is zero. The body motion is reflected in the far field boundary
conditions, and the no-slip boundary condition at the wing is enforced
explicitly through the vorticity and stream function boundary conditions. More
specifically, on the wing, we set Ψ=c, where c is a constant, to satisfy
the no-penetration boundary condition, and ∂Ψ/∂*n*=0 to
satisfy the no-slip boundary condition. At far field, we setΔ×Ψ
=–(**U**_{0}+**r**×Ω_{0}),
where **U**_{0} and Ω_{0} are the translational and
rotational velocity of the wing, respectively, **r** is the position
relative to the wing center, and ω=0. The exact boundary condition onΨ
can be recovered by solving the Poisson equation twice
(Wang, 1999). For this
computation, the far field boundary condition is correct to the dipole
order.

A fourth-order Runge–Kutta scheme is used for the time iterations,
which exhibits a stability domain for this explicit scheme. The stability
condition includes two Courant–Friedrich–Levy (CFL)-like
conditions related to the convection and diffusion time scales over a mesh
size:
3
4
where d*s*=min(dμ, dθ), μ=μ_{0} at the ellipse,
and *C*_{1}=*C*_{2}=0.8. The time step is chosen
to be min(d*t*_{1}, d*t*_{2}). The basic time
iteration in each computation step involves the following sequence:ω
^{n}→Ψ^{n+1}→**u**^{n+1}→ω^{n+1},
where superscripts indicate the time-step. To resolve the flow, 10 grid points
are typically needed along the radial direction in the boundary layer, and at
least 30 points in the azimuthal direction around each tip, whose length scale
is estimated by its radius of curvature. The resolution for
*R*ê100 is 128×256 for the following computations. The
radius of the computational boundary is typically 10 times the chord
length.

The forces on the ellipse can be computed from integrating the stress
tensor along the body. Writing the Navier–Stokes equation in the
coordinates fixed to the wing, we have:
5
6
7
where **U**_{0} and Ω are the translational and rotational
velocity of the wing, and *p* the pressure. The last three terms
corresponds to the non-inertial force due to rotational acceleration, the
Coriolis force, and the centrifugal force. The Coriolis force and the
centrifugal force disappear in the 2D vorticity equation because they can be
recast in terms of the gradient of a potential function.

The velocity and vorticity are solved in the non-inertial coordinates,
which are then transformed into the inertial frame. The forces are calculated
in the inertial frame by integrating the viscous stress:
8
9
where **F**_{p} and **F**_{v} denote pressure and
viscous forces. ρ is the fluid density, *A*_{w} the total
area of the wing, *ŝ* is the tangent vector along the ellipse,
and the integral is over the surface of the ellipse. The second term in
**F**_{p} is similar to the buoyancy force associated with
hydrostatic pressure, i.e. fluids accelerate with–
(d**U**_{0}/d*t*).

### Wing motion, choice of parameters and normalization

The wing follows a sinusoidal flapping and pitching motion. Specifically,
the wing sweeps in the horizontal plane and pitches about its spanwise axis
with a single frequency *f*:
10
11
where *x*(*t*) is the position of the center of the wing, andα
(*t*) is the wing orientation with respect to the
*x*-axis. By definition, the translational and angular velocities are
given by U_{0}(*t*)=d*x*(*t*)/d*t* andΩ
(*t*)=dα(*t*)/d*t*. The parameters include
the stroke amplitude *A*_{0}, the initial angle of attackα
_{0}, the amplitude of pitching angle of attack β, the
frequency *f* and the phase difference φ between
*x*(*t*) and α(*t*). In the experiments, such a
motion is confined to an arc about the wing root, and in the 2D computations,
the motion is along a straight line.

The translational motion of the wing is completely specified by two
dimensionless parameters, the Reynolds number,
*Re*=U_{max}*c*/ν=π*fA*_{0}*c*/ν,
and *A*_{0}/*c*, where U_{max} is the maximum
wing velocity, and *c* the chord. In the subsequent studies, we fix
*f* but vary *A*_{0}/*c* and study its effect on
the flow. For clarity, we will report the value of
*A*_{0}/*c* directly instead of *Re.
A*_{0}/*c* varies from 2.8 to 4.8, with resulting
*Re* from 75 to 115, appropriate for fruitflies. Other parametersα
_{0}, β and *f* are fixed to be π/2, π/4 and
0.25 Hz, respectively.

A main variable of interest in this study is the phase delay between
rotation and translation, φ, which was shown to be a sensitive parameter
in force generation (Dickinson et al.,
1999; Wang,
2000b). Three cases, φ=π/4, 0 and –π/4,
corresponding to the advanced, symmetrical and delayed rotation
(Dickinson et al., 1999), will
be studied for each *A*_{0}/*c*.

We normalize the computational and experimental forces by the maxima of
their corresponding quasi-steady forces, as described in the next section. Our
choice for the normalization is dimensionally the same as the conventional
choice, Gρ*u*^{2}_{rms}*c*, but has the
advantage of normalizing away features specific to the wing, such as its
thickness and geometry. This is because that force dependence of the wing
geometry is sometimes relatively simple. For example, the force coefficients
of ellipses of different thickness were shown to have almost the same
functional dependence on the angle of attack but different magnitude (see
fig. 6 in
Wang, 2000a). The experimental
force coefficients of the robotic fly wing also show little dependence on the
wing planform. If the dependence on the geometry in the steady and unsteady
forces is similar, then their ratio does not depend sensitively on the
geometry of the wing. This would allow us to compare wings of different cross
sections and planforms. Although comparing the force coefficients appears to
be a natural thing to do, one must be cautious when comparing 2D and 3D force
coefficients. The lift coefficient of 1 has different meanings in experiments
and computations, unless the sectional lift coefficient in a 3D wing is a
constant. Strictly speaking, the numbers should only be compared within the
computations or the experiments.

### Quasi-steady forces

Before discussing the unsteady forces, we first describe the calculation of
the quasi-steady forces based on both the translational and rotational
velocity. Because the wing operates at a large range of angles of attack, from
0° to 135°, the Kutta–Joukowski lift, which works for attached
flow associated with small angles of attack, is clearly inapplicable. Instead,
we determine the quasi-steady coefficients empirically, using both the robofly
experiments and computation of a steady translating wing at a fixed angle of
attack. The lift and drag coefficients, defined with respect to the far field
flow, are measured at a time when the forces reach a temporary plateau after
the initial transients (see for example,
fig. 2 in
Dickinson et al., 1999;
fig. 5 in
Wang, 2000a). Forces at all
angles are measured at a fixed time, *t*=2 in dimensionless time
scale.

From the 3D experiments, the lift and drag coefficients are well
approximated by:
12
13
where the angles are expressed in degrees. A similar fit is obtained using our
2D computed data:
14
15
where α is the angle of attack. This fit is shown in
Fig. 1. The constants depend on
the Reynolds number, details of the wing, etc. Unlike the
Kutta–Joukowski lift, which is valid at small values of α and is
proportional to sinα, Equations
12,13
and
14,15
are valid for all values of α and depend explicitly on 2α, rather
than α. The 2α dependence is consistent with the symmetry of a
plate. Moreover, it is consistent with the theoretical prediction of lift and
drag of a stalled wing:
16
17
which is derived assuming complete flow separation in the wake
(von Karman and Burgers,
1963). The theory underpredicts the magnitude of the forces, but
gives roughly the right shape of the force curve. In addition, Equations
16,17
make it apparent that the net force in the stalled case is normal to the wing,
a prediction confirmed by our computations and experiments (Figs
2,
3,
4). We refer to the forces,
Gρ*u*^{2}*C*_{L} and
Gρ*u*^{2}*C*_{D}, as the quasi-steady
translational lift (*L*_{T}) and drag (*D*_{T}),
respectively.

In the case of a translating and rotating wing, the instantaneous velocity
of the wing varies along the chord, as
**u**(**x**,*t*)=**u**_{0}(*t*)+Ω_{0}(*t*)×**x**,
where **x** is the position on the chord measured from the pitching axis.
In the case of constant small pitching amplitude and constant translating
velocity, the potential theory (Munk,
1925) predicts the associated lift to be:
18

Note that Equation 18
includes both the pressure lift of a translating and rotating wing in the
absence of circulation (Magnus force) and the lift due to circulation given by
the Kutta's condition. There is no *a priori* reason as to which of the
quasi-steady forms should fit in the unsteady case better, since both use
assumptions that are invalid in unsteady cases. Recently, revised quasi-steady
models have been proposed to fit these forces
(Sane and Dickinson, 2002).
For the purpose of this paper, we simply compare *L*_{R},
*L*_{T} and *D*_{T} with the unsteady forces. It
turns out that for the prescribed motions here, *L*_{R}
deviates substantially from the unsteady forces, while *L*_{T}
approximates the unsteady forces reasonably well. Therefore, in the subsequent
discussions, we will use *L*_{T} as an estimate for the
quasi-steady forces.

## Results

### Comparison of measured and computed forces

Twelve kinematics are analyzed here: four stroke amplitudes, 60°,
80°, 100° and 120°, at three phase delays, φ=π/4, 0 and–π
/4. The range of amplitudes corresponds to the end-to-end
amplitude:chord ratios of 2.8, 3.6, 4.2 and 4.8 at 70% span. The computational
parameters are chosen to match both the Reynolds number *Re* and the
amplitude:chord ratio, *A*_{0}/*c*. In particular,
*A*_{0} is estimated to be the projection of the arc length at
70% span onto a 2D plane. The location at 70% span was chosen because the
sectional circulation is near maximal there (J. M. Birch, W. Dickson and M. H.
Dickinson, manuscript submitted for publication) and the interference from the
trailing edge vortex is small. Once we fix the computational units for time
and length, i.e. *T*=1 and *c*=2.0, the viscosity is also fixed
in computational units.

Figs 2, 3, 4 summarize the results for variation in φ. Each figure shows the comparison of experimental and computational force coefficients over the first four cycles. In addition, instantaneous force vectors are shown in superposition with the traveling wing during the second stroke.

The wing motion in these cases differs in the angle of attack at the end of stroke. The angles of attack are π/4, π/2 and 3π/4, respectively, in the advanced (φ=π/4), symmetrical (φ=π/2) and delayed (φ=–π/4) rotation cases, as shown in the force vector plots. The delayed rotation case is unusual from the point of view of operating an airfoil. After each reversal, the wing has angles of attack greater thanπ /2, which leads to a downward lift (see Fig. 4). However, insects or bio-mimetic devices may use such a stroke to reduce the force on one wing, and thus generate a torque to turn. In addition, when the wing is moving at an angle greater than π/2, the flow separates quickly, which is qualitatively different from the other two cases. Thus it also provides a good case for testing computations and experiments in different scenarios.

In all cases, the forces quickly settle into an almost periodic state after two strokes. The computational drag follows the experimental drag closely in all three cases. Lift agrees well in the first two cases (Figs 2, 3), but shows a clear phase delay in the case of delayed rotation (Fig. 4). Notice that the shift occurs only after the first stroke. The averaged experimental lift and drag coefficients are (0.93, 1.28), (0.86, 1.34) and (0.38, 1.10), for the advanced, symmetrical and delayed rotation, respectively. The averaged computational lift and drag coefficients are (1.10, 1.36), (0.82, 1.44) and (0.19, 1.21), respectively, for the corresponding cases. We will return to the presence and absence of the phase shift in lift in these three different cases when we discuss the 3D effects.

The averaged force coefficients depend weakly on stroke amplitude, as shown
in Fig. 5. In the case ofφ
=π/4, the average experimental lift coefficients are 0.93, 0.99, 0.95
and 0.93 at *A*_{0}/*c*=2.8, 3.6, 4.2 and 4.8,
respectively. The corresponding computational lift coefficients are 1.07, 1.0,
0.9 and 0.9. The drag coefficients are 1.28, 1.19, 1.12 and 0.93 in
experiments, and 1.36, 1.34, 1.24 and 1.16 in computation. This indicates that
the total force scales roughly with , as
expected. Within this range of amplitude variation, the flows are
qualitatively similar for a given choice of φ. Although we only show theφ
=π/4 case, the results are similar for φ=–π/4 andφ
=0.

The average lift depends sensitively on φ, as emphasized before
(Dickinson et al., 1999;
Wang, 2000b). For example, in
the case of *A*_{0}/*c*=2.8, the averaged values for
experimental lift, *C*_{L} are 0.93, 0.86 and 0.38, for φ
values of 0.25π, 0 and –0.25π, respectively. The comparable
quasi-steady lift coefficients are 0.75, 0.95 and 0.75. This dependence onφ
can be understood intuitively, based on two facts. First, the deviation
between the unsteady forces and quasi-steady forces occurs mostly after the
flip of each stroke. Second, the instantaneous forces in all these cases are
typically normal to the wing, as indicated in the force vector plots in Figs
2,
3,
4, and as discussed in
Materials and methods. Therefore, the contribution to lift and drag can be
correlated with the instantaneous orientation of the wing at the end of each
stroke: π/4, π/2 and 3π/2, for φ=π/4, φ=0 andφ
=–π/4, respectively. One expects an increase in both lift and
drag when φ=π/4, a decrease of lift and increase of drag whenφ
=–π/4, and relatively small change in lift, but a large
increase in drag for φ=0. These indeed are consistent with both the
experimental and computational forces.

### Comparison of experimental and computed vorticity fields

As a further comparison between computation and experiments, we show side by side the snapshots of the vorticity field near the wing from experimental DPIV measurements and computational results (Fig. 6). Ten frames are shown for each period. The colors indicate the strength of the vorticity field. In Fig. 6, columns A and C are computed vorticity, and B and D are experimental vorticity in a 2D slice. The simulations appear to capture the major features of vortex dynamics through a complete stroke cycle. Notice that the fluid momentum is directed downward by pairs of vortices, similar to those shown in asymmetric strokes that model dragonfly wing kinematics (Wang, 2002b). In 3D, the pairs of vortices can be cross-sections of a donut-shaped vortex ring. The structure of the downward momentum jet, characterized by the averaged velocity field over a cycle, is examined elsewhere for both the symmetric and asymmetric strokes strokes (Z. J. Wang, manuscript submitted for publication). Also notice that even the kinematics of left and right strokes are identical, but the flow field differs slightly. This can be seen by comparing columns A and C in Fig. 6. The wing positions are mirror images, but the flows deviate slightly from the mirror symmetry. Such a deviation may be inconsequential in terms of average lift, but worth keeping in mind when interpreting the precise time course of forces.

*Unsteady forces* vs *quasi-steady forces*

The differences between the unsteady and quasi-steady forces have been analyzed extensively for fruitfly kinematics, based on results from tethered animals, with relatively constant translational velocity in the mid-stroke and large acceleration and sharp rotations at the end of strokes. The unsteady effects were dominant near the wing reversal, where they contribute to the rotational effect and the wake capture (Dickinson et al., 1999). The discrepancy between experimental measures of forces and flows (Birch and Dickinson, 2003) and a CFD model of nearly identical conditions (Sun and Tang, 2002) raises a debate about the physical basis of these unsteady effects. Here we do not attempt to resolve these discrepancies, but probe the presence of unsteady effects in a different set of kinematics. In Fig. 7, we compare the unsteady forces with the steady state forces based on the translational velocity. In the case of advanced rotation, the unsteady effects can contribute an additional 50% to the total lift, and in the case of delayed rotation, they can reduce the total lift by a factor of 2–3. It is also clear that the quasi-steady forces based on translational velocity alone do not predict the time-dependent forces during a sinusoidal flapping.

To gauge the relative importance of different unsteady effects, we examine the difference between unsteady forces and estimates based on the translational velocity, as shown in Fig. 7D–F. Ideally we wish to decompose, if possible, the unsteady force approximately into a sum of separate terms, which can be related to wing acceleration, the coupling between rotation and translation, wing–wake interaction, etc. However, in the absence of quantitative prediction of these effects, we can only offer a plausible decomposition by correlating the force peaks with the time course of translational and rotational velocity.

The coupling between translation and rotation can be modeled by
*C*_{rot}Ω(*t*)*u*(*t*), a form
predicted by classical theory of a translating and pitching motion
(Munk, 1925), and tested in
robotic fruitfly experiments (Sane and
Dickinson, 2002), where *C*_{rot} is assumed to be
a constant that depends on the center of rotation. Part of
*C*_{rot}Ω(*t*)*u*(*t*) is the
Magnus force caused by the pressure difference due to velocity difference,
given by Bernoulli's law, and another part is due to additional circulation
caused by the rotational motion to satisfy the Kutta condition
(Munk, 1925). In the three
kinematic patterns studied here, the peaks (labeled `r' in
Fig. 7) associated with
rotation, are picked by matching (with a small shift) the force curve to the
maximum of Ω(*t*)*u*(*t*). The positions of these
measured force peaks vary in the three cases, in accordance with the shift of
the peak positions of *u*(*t*)Ω(*t*). The variation
occurs at the same time scale as
*u*(*t*)Ω(*t*).

Other unsteady effects occur near the wing reversal (labeled `u' in
Fig. 7). The position of these
effects occur roughly at the same time in all three kinematics. These force
peaks do not follow the trace of d*u*(*t*)/d*t*, thus do
not behave as the classical added mass. These force peaks are likely related
to the unsteady growth of vorticity and wake–wake interaction, which do
not have simple analytic expressions in general. Regardless of its physical
basis, the most substantial contribution of this unsteady effect is on drag
(Fig. 7B,C,E,F).

We also note that the peaks alternate in size from stroke to stroke in the
experimental lift, most obvious in Fig.
7A,D. One possible explanation is that this asymmetry reflects the
mechanical artifact due to gear backlash. However a small degree of asymmetry
is also observed in the computational data, e.g. the vorticity field as shown
in Fig. 6 and the forces in
Figs 2,
3,
4. At *Re*=0, the left
and right strokes, which are mirror images about the vertical axis, would
generate forces that have the same symmetry; i.e. the lift from left and right
strokes are identical and the drag forces are of equal size but in the
opposite direction. Here, the Reynolds number is finite and sufficiently large
that the force can depend on the history of flow. One possibility for breaking
the symmetry in forces is by the initial condition.

While identification of the above features is useful when dissecting the
unsteady force traces, it is also relevant to determining their net
contributions. The force directly proportional to the linear acceleration can
have sharp peaks, but it has a zero contribution in reciprocal motions. The
pressure force of a rotating and translating plate, approximated by
*u*(*t*)Ω(*t*), has a non-zero net contribution in
the cases where φ↑0, since
[sin(2π*ft*)cos(2π*ft*+φ)]∼sinφ. The
unsteady vortex force due to wing acceleration has eluded simple analytical
expressions, except for power-law start up flow, where both the added mass
term and the vortex force are calculated analytically and numerically
(Pullin and Wang, 2003). The
unsteady forces contribute to both lift and drag, both predicted in theory
(Pullin and Wang, 2003) and
seen here in Fig. 7.

### Fluid forces and wing inertia

Among various terms contributing to the fluid forces, the pressure force
dominates. The viscous force is smaller by roughly a factor proportional to
1/√*Re*. The pressure force due to non-inertial effects resulting
from translational and rotational acceleration averages to zero in hovering
when the pitching axis is centered at the chord. The magnitude of the
instantaneous non-inertial translational force is also small for these
sinusoidal motions. In particular:
19
where *F*_{ω} is the pressure force due to vorticity,
*F*_{B} is the effective buoyancy due to wing acceleration,
*b* is the thickness of the wing and *A*_{0} is the
amplitude of the stroke. In the derivation of Equation
19, we have used the fact that
dU_{0}/d*t*∼2π*f*U_{0} and
U_{0}∼2π*fA*_{0}.

Fig. 8 illustrates the
contribution of *F*_{ω} (broken line) to the total force
(solid line). In kinematics where there is a fast acceleration at the end of
the stroke, the force will have sharp peaks at the end of the stroke, but the
net contribution is zero, as discussed before.

Finally, we estimate the inertial force associated with wing acceleration
with respect to the fluid force. The inertial force in the experiment turns
out to be negligible compared to the fluid forces, as shown experimentally
(Sane and Dickinson, 2001).
Here is a simple estimate to explain why this is so:
20
where *C*_{f} is the force coefficient, of order 1, *b*
is the thickness of the plate and *A*_{0} the amplitude of the
stroke. For the robofly wing,
(ρ_{wing}/ρ_{fluid})∼1.3 and
(*b*/*A*_{0})≪1, hence the wing inertia is
negligible. For a real insect wing,
(ρ_{wing}/ρ_{fluid})∼10^{3}, so that even
though the wing is very thin, its inertia may not be negligible. Weis-Fogh's
early data showed that the ratio is about 30%
(Weis-Fogh, 1973). It will be
interesting to check this against the biological data. Similar results can be
obtained for the pitching torque:
21
where *c* is the chord. Again, in the experiments the moment of wing
inertia is also negligible during a sinusoidal motion as studied here.

## Discussion

Given that flow around real and model insect wings exhibits 3D effects such as spanwise flow, it is tempting to conclude that 2D computations have little to offer. Here we see that the success and failure of a 2D model in capturing the forces in 3D experiments can provide important insights. In both the advanced and symmetrical rotation cases, the 2D forces are very similar to the 3D forces. Why do they agree? In addition, what do we learn from comparisons between the 2D computation and the 3D experiments?

First, a notable difference between the experimental and computational forces is seen in the delayed rotation, where there is a clear phase shift between the computed and measured lift. In the canonical example of flow past a 2D cylinder and a 3D sphere, the forces during von Karman vortex shedding also show a phase shift, which was argued to be a generic feature between 2D and 3D flows (Mittal and Balachandar, 1995). In view of this, the absence of the phase shift in the advanced and symmetrical cases are particularly interesting. The difference in flow structure in the three cases may be worth further investigation.

Second, these results are relevant to recent discussions about the role of
3D effects on delayed stall. Insects are known to flap their wings at angles
of attack much higher, around 35°
(Ellington, 1984), than the
stall angle of a conventional airfoil, about 12°. As suggested by
Ellington et al. (1996), the
pressure gradient from root to tip within the vortex core might drive spanwise
flow that stabilizes the leading edge vortex by convecting away the vorticity.
The spanwise flow was indeed seen by smoke visualization in the robotic
hawkmoth experiment, where *Re*≈5000
(Ellington et al., 1996;
Willmott et al., 1997). This
proposed mechanism is thought to be analogous to that occurring on delta wing
aircraft, in which spanwise flow through the vortex core maintains the
stability. But as discussed previously
(Ellington et al., 1996), the
exact conditions for establishing spanwise flow in leading-edge vortex for
rotary wings are not completely understood. For example, a helicopter rotor
also experiences a pressure gradient, centrifugal and Coriolis forces, but no
large-scale spanwise flow is observed
(Harris, 1966). Recent smoke
visualization of free-flying butterflies also did not observe substantial
spanwise flow, but reported high variability of 3D flow patterns
(Srygley and Thomas, 2003).
DPIV images of flow field in a robotic fruitfly experiment, where
*Re*≈150, showed no substantial spanwise flow inside the core of
leading edge vortex, but instead indicted substantial spanwise flow behind the
leading edge vortex, which connects to the tip vortex
(Birch and Dickinson, 2001).
Strictly speaking, there is no contradiction among these experiments regarding
the spanwise flow. It is likely that the spanwise flow within the vortex core
occurs only at sufficiently large Reynolds number as in the case of hawkmoth,
but not at low Reynolds number, as in the case of fruitflies. The details of
the spanwise flow can also depend on the wing shape and kinematics. The high
variability of the 3D flow patterns shown by these different experiments,
however, makes it difficult to conclude that spanwise flow is crucial for
generating sufficient lift by a hovering insect.

An alternative explanation for why the conventional stall does not seem to
affect a flapping insect wing relates to the time scale governing the flow
separation that leads to stall. For example, a 2D translating wing at an angle
of attack of 40° and *Re*=1000, does not show a drop in
time-dependent force until the chord travels for about 4 chords, after which
the forces become oscillatory due a von Karman shedding
(fig.∼6 in
Wang, 2000a). Therefore, in
theory there is no need for additional mechanisms to stabilize the leading
edge vortex if the wing travels less than about 4 chords. The early data
compiled by Weis-Fogh (1973)
showed that ratios of stroke-arc to wing-chord of different species during
hovering, including bats (*Plecotus auritus*), birds (hummingbirds),
butterflies and moths (Lepidoptera), wasp and bees (Hymenoptera) and flies
(Diptera), have values less than 4. The beetles (Coleoptera) have values
between 5 and 6. A main difference between a 3D revolving wing and a 2D
translating wing, as noted in recent literature, is that a revolving wing does
not appear to shed its leading edge vortex after a distance much longer than
the stroke length of a typical insect
(Dickinson et al., 1999;
Usherwood and Ellington,
2002). No doubt such a difference would affect the force
coefficients observed in 2D and 3D in the steady state. On the other hand, the
difference in terms of vortex shedding may not affect the transient values. It
is worth re-examining the results of 3D experiments on a flapping wing
(fig. 2D in
Dickinson et al., 1999), which
show that while the 2D steady state lift coefficients underpredict
substantially their the 3D counterparts, the 2D transient values follow
closely the 3D coefficients, up to an angle of attack of about 72°. The 3D
steady force is slightly lower than the unsteady 2D counterpart, due the
well-known downwash due to tip vortices. Similarly, in the cases studied here,
the chord is moving between 3 to 5 chords, and the leading edge vortex does
not appear to separate during each stroke in the cases of advanced and
symmetrical rotation, as indicated by the absence of phase shift between the
2D and 3D forces. In these cases, the 2D forces are good approximations of 3D
experiments.

## ACKNOWLEDGEMENTS

Z.J.W. acknowledges support from AFOSR, NSF Early Career Award and the ONR Young Investigator Program. M.H.D. thanks the Packard Foundation and the National Science Foundation (IBN-0217229) for support of this work.

- © The Company of Biologists Limited 2004