Welcome to our new website

Unsteady forces and flows in low Reynolds number hovering flight: two-dimensional computations vs robotic wing experiments
Z. Jane Wang, James M. Birch, Michael H. Dickinson


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.


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 m3 of mineral oil (density 0.88×103 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 cm2 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.65R, 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.65R as our point of measurement because in a prior DPIV study in which the wake was viewed from the rear, 0.65R 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: Math1 Math2 where u is the velocity field, ω the vorticity field, v is the kinematic viscosity and S is the local scaling factor, S(μ,θ)= cosh2μ–cos2θ 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Δ×Ψ =–(U0+r×Ω0), where U0 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: Math3 Math4 where ds=min(dμ, dθ), μ=μ0 at the ellipse, and C1=C2=0.8. The time step is chosen to be min(dt1, dt2). The basic time iteration in each computation step involves the following sequence:ω n→Ψn+1un+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: Math5 Math6 Math7 where U0 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: Math8 Math9 where Fp and Fv denote pressure and viscous forces. ρ is the fluid density, Aw 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 Fp is similar to the buoyancy force associated with hydrostatic pressure, i.e. fluids accelerate with– (dU0/dt).

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: Math10 Math11 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 U0(t)=dx(t)/dt andΩ (t)=dα(t)/dt. The parameters include the stroke amplitude A0, 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=Umaxc/ν=πfA0c/ν, and A0/c, where Umax is the maximum wing velocity, and c the chord. In the subsequent studies, we fix f but vary A0/c and study its effect on the flow. For clarity, we will report the value of A0/c directly instead of Re. A0/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 A0/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ρu2rmsc, 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.

Fig. 6.

Vorticity plot in the case of A0/c=4.8,φ =0. Ten frames are shown in the fourth stroke. Red, counterclockwise rotating vortices; blue, clockwise rotating vortices. The wing is in black. (A,C) Computed vorticity; (B,D) digital particle image vorticity data in a 2D slice at 0.65R. See Materials and methods for details. Each pair in A,B and C,D corresponds to the same time during a stroke. The time sequence is indicated by the numbers on each plate. The color scale for vorticity of computation and experiments do not correspond to the exact same contour values, so the figure should be viewed as a qualitative comparison.

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.

Fig. 2.

Computational and experimental lift and drag coefficients during advanced rotation (φ=π/4; A0/c=2.8). (A) Lift (CL) and drag (CD) during the first four complete strokes. Red, experimental measurements; blue, computations. The time is normalized with the flapping period. The force is normalized by the maxima of the corresponding quasi-steady forces. (B) Experimental and (C) computational force vectors superimposed on wing positions, plotted at equal time intervals. The green line represents the wing chord; filled circles, the leading edge; arrows indicate force vectors on the wing. R, right; L, left.

Fig. 5.

Force traces for different amplitudes. (A) Experimental force coefficients for advanced rotation (φ=π/4) and three stroke amplitudes, A0/c=2.8 (red), 3.6 (blue) and 4.2 (green). (B) Computational force coefficients for the same parameters. The time is normalized with the flapping period. The force is normalized by the maxima of the corresponding quasi-steady forces.

From the 3D experiments, the lift and drag coefficients are well approximated by: Math12 Math13 where the angles are expressed in degrees. A similar fit is obtained using our 2D computed data: Math14 Math15 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: Math16 Math17 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ρu2CL and Gρu2CD, as the quasi-steady translational lift (LT) and drag (DT), respectively.

Fig. 1.

Quasi-steady lift CL (circles) and drag CD (crosses) coefficients measured from computation compared to the empirical formulae described by Equations 14,15 (solid and broken lines, respectively).

Fig. 3.

Computational and experimental lift and drag coefficients during symmetrical rotation (φ=0; A0/c=2.8). Details as in Fig. 2.

Fig. 4.

Computational and experimental lift and drag coefficients during delayed rotation (φ=–π/4; A0/c=2.8). Details as in Fig. 2.

In the case of a translating and rotating wing, the instantaneous velocity of the wing varies along the chord, as u(x,t)=u0(t)+Ω0(tx, 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: Math18

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 LR, LT and DT with the unsteady forces. It turns out that for the prescribed motions here, LR deviates substantially from the unsteady forces, while LT approximates the unsteady forces reasonably well. Therefore, in the subsequent discussions, we will use LT as an estimate for the quasi-steady forces.


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, A0/c. In particular, A0 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 A0/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 Math, 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 A0/c=2.8, the averaged values for experimental lift, CL 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.

Fig. 7.

Force comparison between experiment, computation and quasi-steady predictions for (A,D) advanced (φ=π/4), (B,E) symmetric (φ=0) and (C,F) delayed (φ=–π/4) rotations. (A–C) Force traces of four strokes starting from rest. Red, experimental force; blue, computational force; green, quasi-steady estimates using Equations 14, 15, 16, 17. Forces are normalized as in Figs 2, 3, 4. (D–F) Difference trace (red) between the experimental forces and quasi-steady forces. Linear acceleration du(t)/dt (blue), and an estimate of rotational force u(t)Ω(t) (green) are shown in arbitrary scale since we are only interested in their basic time course. Broken vertical lines mark the wing reversal during t=[1,2]. Features associated with rotation (r) and unsteady circulation due to wing acceleration (u), are labelled accordingly.

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 CrotΩ(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 Crot is assumed to be a constant that depends on the center of rotation. Part of CrotΩ(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 du(t)/dt, 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: Math19 where Fω is the pressure force due to vorticity, FB is the effective buoyancy due to wing acceleration, b is the thickness of the wing and A0 is the amplitude of the stroke. In the derivation of Equation 19, we have used the fact that dU0/dt∼2πfU0 and U0∼2πfA0.

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.

Fig. 8.

Contributions to drag coefficients. Solid line, total drag; broken line, contribution due to vorticity flux alone. The parameters in the wing kinematics are: A0/c=2.8, φ=π/4.

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: Math20 where Cf is the force coefficient, of order 1, b is the thickness of the plate and A0 the amplitude of the stroke. For the robofly wing, (ρwingfluid)∼1.3 and (b/A0)≪1, hence the wing inertia is negligible. For a real insect wing, (ρwingfluid)∼103, 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: Math21 where c is the chord. Again, in the experiments the moment of wing inertia is also negligible during a sinusoidal motion as studied here.


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.


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.


View Abstract