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

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

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

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

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

**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

*f*:

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

*u*

^{2}

*C*

_{L}and Gρ

*u*

^{2}

*C*

_{D}, as the quasi-steady translational lift (

*L*

_{T}) and drag (

*D*

_{T}),respectively.

**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:

Note that Equation 18includes 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

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

*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:

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

*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:

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

## References

**Birch, J. M. and Dickinson, M. H.**(

**Birch, J. M. and Dickinson, M. H.**(

**Dickinson, M. H. and Götz, K. G.**(

**Dickinson, M. H., Lehmann, F. O. and Sane, S. P.**(

**E, W. and Liu, J.**(

**Ellington, C. P.**(

**Ellington, C. P., van den Berg, C., Willmott, A. P. and Thomas,A. L. R.**(

**Fry, S. N., Sayaman, R. and Dickinson, M. H.**(

*Drosophila.*

**Harris, F.**(

**Liu, H., Ellington, C., Kawachi, K., van den Berg, C. and Willmott, A. P.**(

**Mittal, R. and Balachandar, S.**(

**Munk, M. M.**(

**Pullin, D. I. and Wang, Z. J.**(

**Ramamurti, R. and Sandberg, W. C.**(

**Russell, D. and Wang, Z. J.**(

**Sane, S. and Dickinson, M. H.**(

**Sane, S. and Dickinson, M. H.**(

**Srygley, R. B. and Thomas, A. L.**(

**Sun, M. and Tang, J.**(

**Usherwood, J. R. and Ellington, C. P.**(

**von Karman, T. and Burgers, J. M.**(

**Wang, Z. J.**(

**Wang, Z. J.**(

**Wang, Z. J.**(

**Weis-Fogh, T.**(

**Willmott, A. P., Ellington, C. P. and Thomas, A. L. R.**(

*Manduca sexta.*