Unsteady aerodynamic force generation by a model fruit fly wing in flapping motion

It has been shown that conventional aerodynamic theory, which was based on fixed-wing aircraft and steady-state flow conditions, cannot explain the generation of large lift by the wings of small insects (for reviews, see Ellington, 1984a; Spedding, 1992). In the past few years, a great deal of work has been conducted to reveal the unsteady mechanisms involved. Dickinson and Götz (1993) measured the aerodynamic forces of an aerofoil started impulsively (Reynolds number Re=100) and showed that the lift was enhanced by the presence of a dynamic stall vortex. But lift enhancement was limited to only approximately 2–3 chord lengths of travel because of the shedding of the dynamic stall vortex. The wings of most insects travel approximately five chord lengths during an up- or downstroke (Weis-Fogh, 1973).

Ellington et al. (1996) performed flow-visualization studies on a hawkmoth Manduca sexta using a mechanical model of the hawkmoth wings. They discovered that the dynamic stall vortex on the wing did not shed during the translational motion of the wing in both the up- and downstrokes because it was stabilized by a strong spanwise flow. The authors suggested that this was a new mechanism of lift enhancement, prolonging the benefit of the delayed stall for the entire stroke. This mechanism of lift enhancement was confirmed by computational fluid-dynamic analyses (Liu et al., 1998; Lan and Sun, 2001).

Recently, Dickinson et al. (1999) conducted force measurements on flapping robotic fruit fly wings and showed that a large lift force was maintained during the translational phases of the up- and downstrokes and brief bursts of high lift occurred during stroke reversal when the wing was rotating. They explained the large lift force during the translational phases by the delayed-stall mechanism of Ellington et al. (1996) and suggested some possible explanations for the large lift force recorded during wing rotation. Only force measurements were conducted by Dickinson et al. (1999); if simultaneous flow-field information had been obtained, further insights into the high lift generation process could be obtained.

In the present paper, we conduct a computational fluid-dynamic analysis on unsteady aerodynamic force generation for a model fruit fly wing in flapping motion. Unsteady aerodynamic forces and flow fields are obtained simultaneously by numerical solution of the Navier–Stokes equations. The high lift generation process can be explained further on the basis of force and flow-field information. The flapping motion of normal hovering flight (Fig. 1), similar to that used by Dickinson et al. (1999), is considered.

As in Dickinson et al. (1999), the planform of the wing considered in the present study is approximately the same as that of a fruit fly (Drosophila sp.) wing (see Fig. 2B). The wing section is an ellipse whose thickness is 12 % of the aerofoil chord length, and the radius of the leading and trailing edges is 0.2 % of the aerofoil chord length. Following biomechanics convention, the azimuthal rotation of the wing about the y axis is called ‘translation’ and the pitching rotation of the wing near the end of a stroke and at the beginning of the following stroke is called ‘rotation’. The speed at the span location r₀ is called the translational speed of the wing, and r₀ is the radius of the second moment of wing area and is determined by r₀=(∫_Sr²dS/S)^½, where r is the radial distance and S is the wing area. For the model fruit fly wing considered, the distance between the wing root and the y axis is 0.36c, where c is the mean chord length, the wing-span is 2.4c, the total wing length l (the distance between the y axis and the wing tip) (Fig. 1A) is 2.76c, and r₀ is 1.6c or 0.58l.

The flapping motion considered in the paper is an idealized one, similar to that considered by Dickinson et al. (1999). A stroke consists of the following three parts, as shown in Fig. 1B: pitching-down rotation and translational acceleration at the beginning of the stroke, translation at constant speed and constant angle of attack during the middle of the stroke, and pitching-up rotation and translational deceleration at the end of the stroke. In normal hovering flight, the wing motion during the upstroke is identical to that during the downstroke. The translational speed is denoted by u_t, which takes a constant value U (the reference velocity) except at the beginning and near the end of a stroke. During the deceleration at the end of a stroke and the acceleration at the beginning of the following stroke, u_t is given by:

where u_t⁺=u_t/U, τ=t*U/c, where t* is (dimensional) time, τ₁ is the time at which the deceleration near the end of a stroke starts and τ₁+Δτ_t is the time at which the acceleration at the beginning of the next stroke finishes. In the calculation, u_t⁺ determines the azimuthal-rotational speed of the wing. Denoting the azimuthal-rotational speed as ψ̇, we have ψ̇(τ)=u_t/r₀, where τ is non-dimensional time. The angle of attack of the wing is denoted by α, which also takes a constant value except at the beginning or near the end of a stroke. Following the experimental work of Dickinson et al. (1999), α is set as 40° in the present study. Around stroke reversal, α changes with time, and the angular velocity, α̇, is given by:

where α̇⁺= α̇c/U, α̇₀⁺ is a constant, τ_r is the non-dimensional time at which the rotation starts and Δτ_r is the time interval over which the rotation lasts. In the time interval Δτ_r, the wing rotates from α=40° to α=140°. Therefore, when Δτ_r is specified, α̇₀⁺ can be determined (around the next stroke reversal, the wing would rotate from α=140° to α=40°, so the sign of the right-hand side of equation 2 should be reversed). The Reynolds number is defined as Re=cU/ν, where ν is the kinematic viscosity, and Re=136 is used here, as in Dickinson et al. (1999), as typical of a fruit fly wing. A stroke angle of 155° is assumed, approximately the same as that used by Dickinson et al. (1999). On the basis of the stroke angle, the non-dimensional time taken for one cycle (two strokes) is approximately 10.8. In the calculation, the wing starts the flapping motion in still air. The calculation is stopped when periodicity in the aerodynamic force and flow structures is approximately reached.

In the flapping motion considered in the present paper, the wing performs translational motion (azimuthal rotation) and pitching rotation. To calculate the flow around a body performing unsteady motion (such as the present flapping wing), one approach is to write and solve the governing equations in a body-fixed, non-inertial frame of reference with inertial force terms added to the equations. An advantage of this approach is that the coordinate transformation that generates a body-conforming computational grid does not need to vary with time, and the grid is generated only once. But, in this approach, some treatment is needed to handle the far-field boundary conditions when the body is rotating (velocity at the far-field boundary tends to be infinite), introducing extra terms into the equations.

Another approach is to write and solve the governing equations in an inertial frame of reference. By using a time-dependent coordinate transformation and the relationship between the inertial and non-inertial frames of reference, a body-conforming computational grid in the inertial frame of reference (which varies with time) can be obtained from a body-conforming grid in the body-fixed non-inertial frame, which needs to be generated only once. This approach does not need special treatment for the far-field boundary conditions and, moreover, since no extra terms are introduced into the equations, existing numerical methods can be applied directly to the solutions of the equations. This approach is employed here.

The non-dimensionalized three-dimensional incompressible unsteady Navier–Stokes equations, written in the inertial coordinate system o,x,y,z (see Fig. 1A), are as follows:

where u, v and w are three components of the non-dimensional velocity and p is the non-dimensional pressure. In the non-dimensionalization, U, a constant velocity given above, c, the mean chord length of the wing, and c/U are taken as the reference velocity, length and time, respectively. Equations 3–6 are transformed from the Cartesian coordinate system (x,y,z,t) to the curvilinear coordinate system (ξ,η,ζ,τ) using a general time-dependent coordinate transformation in the form:

The transformed equations written in conservative form are as follows:

where J is the Jacobian of the transformation and:

where the symbol ∇ is the gradient operator, and the velocity gradients and the metrics of the transformation were written as:

Fig. 1A also shows a wing-fixed coordinate system (o′,x′,y′,z′). The inertial coordinates (o,x,y,z) are related to the wing-fixed coordinates (o′,x′,y′,z′) through the following relationship:

where a is the distance between o and o′. Using equation 22, the metrics in the inertial coordinate system, (ξ_x,ξ_y,ξ_z,ξ_t), (η_x,η_y,η_z,η_t), (ζ_x,ζ_y,ζ_z,ζ_t), which are needed in equations 8 and 9, can be calculated from those in the body-fixed, non-inertial coordinate system, (ξ_x′,ξ_y′,ξ_z′), (η_x′,η_y′,η_z′), (ζ_x′,ζ_y′,ζ_z′), which need to be calculated only once. As the wing moves, the coordinate transformation functions vary with (x,y,z,t) such that the grid system moves and always fits the wing. The wing-fixed non-inertial frame of reference (o′,x′,y′,z′) is used in the initial grid generation and also in the description of the calculated results.

Equations 8 and 9 are solved using the algorithm developed by Rogers and Kwak (1990) and Rogers et al. (1991) and are in the same form as that solved by Rogers et al. (1991). The algorithm is based on the method of artificial compressibility, which introduces a pseudotime derivative of pressure into the continuity equation. Time accuracy in the numerical solutions is achieved by subiterating in pseudotime for each physical time step. The algorithm uses a third-order flux-difference splitting technique for the convective terms and the second-order central difference for the viscous terms. The time derivatives in the momentum equation are differenced using a second-order, three-point, backward-difference formula. The algorithm is implicit and has second-order spatial and time accuracy. For details of the algorithm, see Rogers and Kwak (1990) and Rogers et al. (1991).

At the inflow boundary, the velocity components are specified as free-stream conditions while the pressure is extrapolated from the interior. At the outflow boundary, the pressure is set equal to the free-stream static pressure and the velocity is extrapolated from the interior. On the wing surface, impermeable wall and no-slip boundary conditions were applied, and the pressure on the boundary is obtained through the normal component of the momentum equation.

A body-conforming grid was generated by using a Poisson solver based on the work of Hilgenstock (1988). The grid topology used was an O-H grid. A portion of the grid used for the wing is shown in Fig. 2.

A code based on the above numerical method was written by Lan and Sun (2001), and it was verified by the analytical solution of the boundary layer flow on a flat plate and validated by comparing the calculated and measured pressure distributions on a wing. When the code is used in the present study, before making observations pertaining to the physical aspect of the flow, estimates of the accuracy of the computed solutions must be provided. In the present calculations, numerical uncertainties are mainly associated with time discretization (i.e. time step values), spatial resolution and far-field boundary location. The effects of these variables on the computed solutions will be discussed below together with the calculated results.

Once the Navier–Stokes equations have been numerically solved, the velocity components and pressure at discretized grid points for each time step are available. The aerodynamic force acting on the wing arises from the surface pressure and viscous stress along the wing surface. Integrating the pressure and viscous stress over the wing surface at a time step gives the total force acting on the wing at the corresponding time instant. The lift of the wing, L, is the component of the total force perpendicular to the translational velocity and is positive when it is directed upwards. The drag, D, is the component of the total force parallel to the translational velocity and is positive when directed opposite to the direction of the translational velocity of the downstroke. The lift and drag coefficients, denoted by C_L and C_D, respectively, are defined as follows:

where ρ is the fluid density.

First, we considered a case (henceforth called case A) with the following conditions: Δτ_r=3.48, approximately the same as in Dickinson et al. (1999); Δτ_t=2.6, approximately twice that used by Dickinson et al. (1999) (smaller values will be used in future cases), the wing-rotation axis is at 0.2c from the leading edge of the wing and most of the wing rotation occurs near the end of a stroke with only a very small portion at the beginning of the next stroke (wing rotation precedes stroke reversal). As found by Dickinson et al. (1991), when wing rotation preceded stroke reversal, a large lift force could be produced. Fig. 3 gives the translational (u_t⁺) and rotational (α̇⁺) velocities of the wing and the lift (C_L) and drag (C_D) coefficients versus τ during one cycle (the time constants, τ₁ and τ_r in equations 1 and 2 can be read from Fig. 3A). The contours of the non-dimensional spanwise component of the vorticity ω_z′ are given in Fig. 4 (in the non-dimensionalization, U and c were used as reference velocity and length, respectively) and sectional streamline plots (seen in the o,x′,y′,z′ frame moving with the wing) in Fig. 5. These figures will be used to explain the aerodynamic force behaviour.

In Fig. 3, the results calculated using three different grids are plotted. Grid 1 had dimensions 45×53×35 in the normal direction, around the wing section and in the spanwise direction, respectively, grid 2 was 62×77×49 and grid 3 was 93×109×71. The normal grid-spacings at the wall for the above three grids were 0.003, 0.002 and 0.002, respectively. The outer boundary for all three grids was set at 10 chord lengths from the wing, and the time step used was 0.02. There was only a slight difference between the results calculated using grids 1 and 2 and almost no difference in the results for grids 2 and 3. Calculations were also performed using a larger computational domain. To isolate the effects of domain size, the outer boundary was sited farther from the wing by adding more grid points in the normal direction of grid 3. The results showed that siting the outer boundary more than 10 chord lengths away from the wing had no effect on the results. It was concluded that a grid size of 93×109×71 was appropriate for the present study.

Variations in C_L and C_D over a flapping cycle in the present study were similar to those of Dickinson et al. (1999) (see their Figs 1D and 3A). A comparison of the calculated lift with the measured values will be given below.

As noted above, since the flow conditions in the down- and upstrokes are essentially the same, only one stroke is discussed here. At the beginning of the stroke (τ=16.7–18), C_L increases to peak at approximately 1.2, and variation in C_D is small. During this period, the wing accelerates towards the right and rotates clockwise by only a small angle, i.e. from α=134.4° at τ=16.7 to α=140° at τ=17.5 (see Fig. 4A,B); thus, its motion is mainly translational acceleration. In fast translational acceleration, as shown by Lan and Sun (2001), a peak in C_L could be produced as a result of the generation of strong vorticity layers in a short period, giving a large rate of change of fluid impulse. As can be seen in Fig. 4C, at the end of the acceleration, a new layer of positive vorticity has formed around the leading edge and upper surface of the wing, and a new layer of negative vorticity has formed on the lower surface of the wing extending beyond its trailing edge (starting vortex).

Between τ=18 and τ=19.4, the wing translates (azimuthally rotates) with constant speed at α′=40° (where α′=180°–α and α=140°); a large C_L, approximately 1.3, is maintained. If the wing is not in azimuthal rotation but is moving like an aeroplane wing, as shown by Lan and Sun (2001), stall would accur at approximately two chord lengths of travel after starting, or at τ≈18.7, and C_L would start to decrease at about this time. In the present case (Figs 4C–E, 5A–C), the dynamic stall vortex does not shed, stall is absent and a large C_L can therefore be maintained; this is the ‘stall-absence mechanism’ revealed by Ellington et al. (1996).

Near the end of the stroke, between τ=19.4 and τ=22.1, large values of C_L and C_D occur, peaking at τ≈20.6 (Fig. 3). From τ=19.4 to τ=20.6, while still translating with constant speed, the wing rotates (pitches up) rapidly. This motion, fast pitching-up while translating with constant speed, results in sharp increases in C_L and C_D, also due to the generation of strong vorticity layers over a short period (Lan and Sun, 2001). As seen in Fig. 4F, a new layer of positive vorticity is formed around the leading edge of the wing and a new layer of negative vorticity near the trailing edge of the wing. From τ=20.8 to τ=22.1, the wing is in deceleration, causing sharp decreases in C_L and C_D (Lan and Sun, 2001). The above discussion shows that pitching-up rotation causes C_L and C_D to increase rapidly, and deceleration causes them to decrease rapidly, forming the large peaks in C_L and C_D in the last part of the stroke.

When the wing moves like an aeroplane wing under steady-state conditions, at the same Reynolds number and angle of attack (Re=136 and α=40°), the calculated C_L is 0.59. Vogel (1967) measured the forces on a thin plate cut to the shape of a fruit fly wing in a wind tunnel, at Re=200 and α=40°, C_L was approximately 0.6. In fact, from our steady-state calculations and from the measurements of Vogel (1967), the maximum value of C_L at steady-state conditions is approximately 0.6. For reference, we used this maximum C_L value and the translational velocity of the wing, u_t (velocity at r₀), to estimate the quasi-steady lift of the wing, 0.6(0.5ρu_t²S), and therefore the quasi-steady lift coefficient, 0.6(u_t²/U²). The quasi-steady lift coefficient is plotted in Fig. 3B for comparison. For the flapping motion described above, C_L is much larger than the corresponding quasi-steady value for most of the stroke. The mean value of C_L is 1.2, more than twice the quasi-steady value (0.5 in Fig. 3). As seen in Fig. 3, this large mean C_L results from three causes: the C_L peak at the beginning of the stroke, the large C_L in the middle of the stroke and the large C_L peak near the end of the stroke. The above analysis showed that the C_L peak at the beginning of the stroke could be explained by the rapid acceleration of the wing, the large C_L in the middle of the stroke by the stall-absence mechanism and the large C_L peak near the end of the stroke by the rapid pitching-up rotation of the wing while in constant-speed translation. However, Dickinson et al. (1999) gave somewhat different explanations for the C_L peaks at the beginning and near the end of the stroke. These explanations are investigated in greater detail below.

Effects of acceleration at the beginning of a stroke and effects of the wake of the previous strokes

Fig. 6 presents results for when Δτ_t is varied in equation 1 (smaller Δτ_t corresponds to a larger acceleration at the beginning of a stroke); other conditions are the same as in case A. A larger acceleration results in a larger C_L peak at the beginning of the stroke, indicating that the C_L peak at the beginning of a stroke is closely related to the rapid acceleration of the wing. Dickinson et al. (1999) considered that the large C_L peak at the beginning of a stroke could be explained by the mechanism of wake capture. This mechanism was revealed by Dickinson (1994) for a two-dimensional aerofoil in a simplified flapping motion with two opposite strokes, in which the flow generated by the first stroke could increase the effective fluid velocity, and hence the lift, at the start of the next stroke. In a recent study by Sun and Hamdani (2001), it was shown that, for a two-dimensional aerofoil, if the first stroke was very short and the dynamic stall vortex did not shed (or a vortex street similar to the von Karman vortex street did not form), the mechanism of wake capture would not exist. For a three-dimensional wing in flapping motion, Ellington et al. (1996) showed (see also above) that the dynamic stall vortex does not shed and, therefore, that the mechanism of wake capture might not exist.

To investigate the effect of the vorticity wake left from the previous strokes, the wing was started from still air and moved in the same way as in a stroke of the flapping motion. The results are shown in Fig. 7 together with the results from Fig. 6 for comparison. The difference between these two sets of results represents the effect of the wake from the previous strokes. At the beginning of the stroke (and throughout the stroke), the effect of the wake from the previous stroke does not increase the C_L, slightly decreasing it instead. This detrimental effect of the wake from the previous stroke on C_L can be predicted from the vorticity plots shown in Fig. 4. At the start of the stroke (Fig. 4A), the dynamic-stall vortex from the previous stroke (dashed lines in Fig. 4A), which has clockwise vorticity, is under the leading edge of the wing and the trailing edge from the last stroke (solid lines in Fig. 4A), which has counterclockwise vorticity, is below the wing. These vortices are positioned such that they would produce downwash velocity in front of the wing, decreasing its lift. Fig. 8 gives the absolute velocity vectors in the middle section of the wing, corresponding to the vorticity field in Fig. 4A. At this point, wing rotation has almost finished and translation has just started, so there is almost no motion of the wing and the flow velocity shown in Fig. 8 is caused by the wake from the previous stroke. The downwash in front of the wing can be clearly seen. The above results show that the wake from the previous stroke does not increase the lift and that the large C_L peak at the beginning of a stroke is due to the rapid acceleration of the wing.

From Fig. 7, it can also be seen that, when the acceleration is larger (smaller Δτ_t, Fig. 7B), the effect of the wake from previous strokes becomes weaker and that, during the rapid pitching-up rotation of the wing near the end of a stroke, the effect of the wake is very small.

The large C_L peak near the end of a stroke was explained by Dickinson et al. (1999) as analogous to the Magnus effect. When a moving cylinder or sphere spins, the friction between the fluid and the surface of the body tends to drag the fluid near the surface in the same direction as the rotational motion. Superimposed onto the usual non-spinning flow, this ‘extra’ velocity contribution creates a higher-than-usual velocity (hence lower pressure) at the top of the body and a lower-than-usual velocity (hence higher pressure) at the bottom, resulting in lift. This phenomenon is called the Magnus effect. The Magnus effect is a steady-state flow phenomenon. Note that it cannot explain the large C_D peak that accompanies the large C_L peak (see Figs 3 and 6) (see also fig. 1D in Dickinson et al., 1999). It was shown above that the large C_L and C_D peaks at the end of the stroke were due to the effects of rapid pitching-up rotation. This effect would be expected to become weaker when the location of the rotation axis is moved rearward or when the rotation rate is lower. Fig. 9 gives the results for three different locations of the rotation axis (including that in case A, where it was located at 0.2c from the leading edge of the wing); other conditions are the same as in case A. When the rotation axis is moved rearward, the peak C_L (and C_D) decreases greatly.

Fig. 10 presents the results for different rotation rates (varying Δτ_r in equation 2 while the other conditions are kept the same) (see Fig. 10A for the motion conditions). A lower rotation rate (larger Δτ_r) produces a smaller C_L peak. These results further show that the large C_L peak near the end of a stroke is due to the pitching-up rotation of the wing.

In the above analyses, wing rotation preceded stroke reversal. Here, the effects of shifting the rotation in time are studied. Two cases are considered. In the first, rotation occurs symmetrically with respect to stroke reversal, i.e. half the rotation is conducted in the latter part a stroke and the other half in the early part of the following stroke. In the second case, rotation is delayed with respect to stroke reversal. Other conditions are as in case A. The lift and drag coefficients are shown in Fig. 11 together with results for case A (rotation advanced) for comparison. The timing of the rotations is given in Fig. 11A.

In the case of symmetrical rotation, compared with the rotation-advanced case (case A), near the end of the stroke the increases in C_L and C_D start later and lower peak values are attained. This is because there is less time for the wing to conduct pitching-up rotation while translating at constant speed. At the beginning of the stroke, C_L is also smaller (no C_L peak exists) because, at this point, the wing is conducting pitching-down rotation, cancelling some of the effect of the translational acceleration. The mean C_L is 0.91 which, although smaller than that for case A, is still 80 % larger than the quasi-steady value.

Fig. 12 gives examples of the vorticity plots for symmetrical rotation. At the beginning of the stroke (Fig. 12A,B), the wing performs the same translational acceleration as in case A (from τ=16.7 to τ=18, towards the right). But since it also rotates clockwise (from α=86.5° to α=140°) and the rotation axis is near the leading edge of the wing (0.2c from the leading edge), a large part of the wing moves backwards (towards the left). As a result, the starting vortex (Fig. 12B) has less vorticity and moves less far downstream from the wing than in case A at the corresponding time (Fig. 4C). Sun and Hamdani (2001) and Lan and Sun (2001) showed that, if the starting vortex is weaker and moves less far downstream in a given time, less lift will be generated. This helps to explain the smaller C_L at the beginning of the stroke. Near the end of the stroke, e.g. at τ=20.9 (Fig. 12D), the wing has just started pitching-up rotation and a small new vorticity is generated [at the corresponding time in case A (see Fig. 4G) a strong new vorticity layer is generated]. This helps to explain the smaller C_L peak near the end of the stroke.

In the rotation-delayed case (Fig. 11), near the end of a stroke, there is no C_L peak; instead, C_L decreases to very small, even negative, values. This is because, in this period, the wing is decelerating rapidly and the pitching-up rotation rate is still very small (the majority of the rotation is delayed until the next stroke). At the beginning of the stroke, although it is in fast acceleration, the wing has a simultaneous rapid pitching-down rotation, resulting in a very small, even negative, C_L. The mean C_L is only 0.50, much less than in case A. Fig. 13 gives examples of vorticity plots for the rotation-delayed case. Note that, at the start of the stroke, the wing has rotated by only a few degrees (Fig. 13A) and that the wing attitude and the vorticity around the wing are very different from those of case A (Fig. 4A). During the translational acceleration at the beginning of the stroke (τ=16.7 to τ=18), the wing rotates clockwise, from α=47.4° to α=115.5° (Fig. 13A,B). The formation of the starting vortex is suppressed by this rotation (compare Fig. 13B with Figs 4C and 12B), and a vortex with positive vorticity (opposite in sign to that of the starting vortex) is formed at the leading edge of the wing. This helps to explain the negative C_L at the beginning of the stroke. In the last part of the stroke, from at τ=20.8 to τ=22.1, the wing is in rapid deceleration and has almost no rotation. As a result, a new layer of positive vorticity on the lower surface and around the trailing edge of the wing and a new layer of negative vorticity on the upper surface and around the leading edge of the wing are formed (compare Fig. 13F with Fig. 13E). Lan and Sun (2001) showed that such formation of new vorticity layers, opposite to that of a wing in acceleration, would cause the lift to decrease sharply. This helps to explain the C_L pattern in the last part of the stroke.

We have presented above a detailed analysis of the unsteady forces that occur in connection with the flow structure. Here, we compare the results of those calculations with experimental results and with data applicable to free-hovering flight of the fruit fly Drosophila virilis.

Dickinson et al. (1999) measured the unsteady lift of a robotic wing modelled on the fruit fly. Their lift is presented in dimensional form and, for comparison, we need to convert it into the lift coefficient. In their experiment, the fluid density ρ was 0.88×10³ kg m^–3, the wing length l was 0.25 m and the wing area S was 0.0167 m². The translational speed at the wing tip during the constant-speed translation phase of a stroke, given in fig. 3D of Dickinson et al. (1999), is 0.235 m s^–1 and, therefore, the reference speed (speed at r₀=0.58l) is calculated as U=0.14 m s^–1. From the above data,

Data for free hovering flight of the fruit fly Drosophila virilis were taken from Weis-Fogh (1973); insect weight was 1.96×10^–5 N, wing length l was 0.3 cm, the area of both wings S_t was 0.058 cm², stroke angle φ was 2.62 rad and stroke frequency n was 240 s^–1.

We investigated whether the insect weight can be balanced by the mean lift calculated in the present study. From the definition of C_L in equation 23:

where L̅ is the mean lift and C̅_L̅ is the mean lift coefficient (air density ρ is taken as 1.25×10^–6 g cm^–4 s²). The speed at the radial location r₀ was defined as the translational speed of the wing, and the translational speed during the constant-speed translation phase of a stroke was taken as the reference speed U. To calculate of the results shown in Figs 3 and 11, the mean translational speed over a stroke was 0.825U. For the hovering fruit fly of Weis-Fogh (1973), the mean translational speed is 2φnr₀=219 cm s^–1, and U can be calculated as U=219/0.825=265 cm s^–1. Inserting the values of ρ, U and S_t into equation 25, the mean lift can be written as

For the symmetrical rotation and rotation-advanced cases, C̅_L̅ was 0.91 and 1.20, respectively. Both exceed that needed to balance the weight. It should be noted, however, that the angle of attack at midstroke used in these calculations was 40°, which may be larger than the actual value. Calculations were performed with smaller angles of attack at midstroke (other conditions kept the same). It was shown that, when the angle of attack was 36°, C̅_L̅ was 0.79 for symmetrical rotation and 1.02 for the rotation-advanced case; when the angle of attack was 34°, C̅_L̅ was 0.75 and 0.98, respectively.

The above results show that, with symmetrical rotation and an angle of attack at midstroke of approximately 36°, the fruit fly can produce enough lift for hovering flight. Ellington (1984b) observed many small insects in hovering flight, including the fruit fly, and found an angle of attack of approximately 35°, so our predicted value is in very good agreement with the observations. The above results also show that, if the insect employs a larger angle of attack or changes the timing of wing rotation, much greater lift can be produced for manoeuvring or other purposes.

a    distance between the origins of the inertial and non-inertial frames of reference

c    mean chord length

C_D    drag coefficient

C_L    lift coefficient

D    drag

J    Jacobian of the transformation between (x,y,z,t) and (ξ,η,ζ,τ)

l    wing length

L    lift

L̅    mean lift

n    wingbeat frequency

o    origin of the inertial frame of reference

o′    origin of the non-inertial frame of reference

p    non-dimensional fluid static pressure

r    radial position along the wing length

r₀    radius of the second moment of wing area

Re    Reynolds number

S    area of one wing

S_t    area of a wing pair

t*    time

u,v,w    non-dimensional velocity component in the x,y,z directions, respectively

u_t    translational velocity of the wing

u_t⁺    non-dimensional translational velocity of the wing

U    reference velocity

x,y,z    coordinates in the inertial frame of reference

x′,y′,z′    coordinates in the non-inertial frame of reference

α    angle of attack

α′    180°–α

α̇    angular velocity of pitching rotation

α̇⁺    non-dimensional angular velocity of pitching rotation

Δτ_r    duration of pitching rotation, non-dimensional

Δτ_t    duration of translational acceleration, non-dimensional

ν    kinematic viscosity

ξ,η,ζ    transformed coordinates

ρ    density of fluid

τ    non-dimensional time (τ=t)

τ₁    time when translational deceleration starts, non-dimensional

τ_r    time when pitching rotation starts, non-dimensional

φ    stroke amplitude

ψ    azimuthal rotation angle

ψ̇    angular velocity of azimuthal rotation

ω_z′    spanwise component of vorticity, non-dimensional

∇    gradient operator (∂/∂x,∂/∂y,∂/∂x)

This research was supported by the National Natural Science Foundation of China. We thank the two referees whose thoughtful questions and valuable suggestions greatly improved the quality of the paper.