A computational fluid-dynamic analysis was conducted to study the unsteady aerodynamics of a model fruit fly wing. The wing performs an idealized flapping motion that emulates the wing motion of a fruit fly in normal hovering flight. The Navier–Stokes equations are solved numerically. The solution provides the flow and pressure fields, from which the aerodynamic forces and vorticity wake structure are obtained. Insights into the unsteady aerodynamic force generation process are gained from the force and flow-structure information.

Considerable lift can be produced when the majority of the wing rotation is conducted near the end of a stroke or wing rotation precedes stroke reversal (rotation advanced), and the mean lift coefficient can be more than twice the quasi-steady value. Three mechanisms are responsible for the large lift: the rapid acceleration of the wing at the beginning of a stroke, the absence of stall during the stroke and the fast pitching-up rotation of the wing near the end of the stroke.

When half the wing rotation is conducted near the end of a stroke and half at the beginning of the next stroke (symmetrical rotation), the lift at the beginning and near the end of a stroke becomes smaller because the effects of the first and third mechanisms above are reduced. The mean lift coefficient is smaller than that of the rotation-advanced case, but is still 80 % larger than the quasi-steady value. When the majority of the rotation is delayed until the beginning of the next stroke (rotation delayed), the lift at the beginning and near the end of a stroke becomes very small or even negative because the effect of the first mechanism above is cancelled and the third mechanism does not apply in this case. The mean lift coefficient is much smaller than in the other two cases.

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 Gö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.

### The model wing and its kinematics

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 r0 is called the translational speed of the wing, and r0 is the radius of the second moment of wing area and is determined by r0=(∫Sr2dS/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 r0 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 ut, 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, ut is given by:

$\mathit{u}_{t}^{{+}}{\ }{=}{\ }0.5{[}1{\ }{+}{\ }cos{\pi}({\tau}{\ }{-}{\ }{\tau}_{1})/{\Delta}{\tau}_{t}{]};\ {\tau}_{1}{\leqslant}{\tau}{\leqslant}{\tau}_{1}{+}{\Delta}{\tau}_{t}\ ,$
1

where ut+=ut/U, τ=t*U/c, where t* is (dimensional) time, τ1 is the time at which the deceleration near the end of a stroke starts and τ1+Δτt is the time at which the acceleration at the beginning of the next stroke finishes. In the calculation, ut+ determines the azimuthal-rotational speed of the wing. Denoting the azimuthal-rotational speed as ψ̇, we have ψ̇(τ)=ut/r0, 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:

${\dot{{\alpha}}}^{{+}}{\ }{=}{\ }0.5{\dot{{\alpha}}}_{0}^{{+}}{\{}1{\ }{-}{\ }cos{[}2{\pi}({\tau}{\ }{-}{\ }{\tau}_{r})/{\Delta}{\tau}_{r}{]}{\}};\ {\tau}_{r}{\leqslant}{\tau}{\leqslant}{\tau}_{r}{+}{\Delta}{\tau}_{r}\ ,$
2

where α̇+= α̇c/U, α̇0+ 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, α̇0+ 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.

### The Navier–Stokes equations and the computational method

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:

3
4
5
6

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

${\xi}{\ }{=}{\ }{\xi}(\mathit{x},\mathit{y},\mathit{z},\mathit{t}){\ },\ {\eta}{\ }{=}{\ }{\eta}(\mathit{x},\mathit{y},\mathit{z},\mathit{t}){\ },\ {\zeta}{\ }{=}{\ }{\zeta}(\mathit{x},\mathit{y},\mathit{z},\mathit{t})\ and\ {\tau}{\ }{=}{\ }\mathit{t}\ .$
7

The transformed equations written in conservative form are as follows:

8
9

where J is the Jacobian of the transformation and:

$\mathit{A}\ {=}\ {\xi}_{x}\mathit{u}\ {+}\ {\xi}_{y}\mathit{v}\ {+}\ {\xi}_{z}\mathit{w}\ ,$
10
$\mathit{B}\ {=}\ {\eta}_{x}\mathit{u}\ {+}\ {\eta}_{y}\mathit{v}\ {+}\ {\eta}_{z}\mathit{w}{\ },$
11
$\mathit{C}\ {=}\ {\zeta}_{x}\mathit{u}\ {+}\ {\zeta}_{y}\mathit{v}\ {+}\ {\zeta}_{z}\mathit{w}{\ },$
12
13
14
15
16
17
18
19

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

20
21

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:

22

where a is the distance between o and o′. Using equation 22, the metrics in the inertial coordinate system, (ξxyzt), (ηxyzt), (ζxyzt), 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.

### Evaluation of the aerodynamic forces

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 CL and CD, respectively, are defined as follows:

23
24

where ρ is the fluid density.

### A typical case

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 (ut+) and rotational (α̇+) velocities of the wing and the lift (CL) and drag (CD) coefficients versus τ during one cycle (the time constants, τ1 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 CL and CD 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), CL increases to peak at approximately 1.2, and variation in CD 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 CL 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 CL, 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 CL 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 CL 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 CL and CD 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 CL and CD, 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 CL and CD (Lan and Sun, 2001). The above discussion shows that pitching-up rotation causes CL and CD to increase rapidly, and deceleration causes them to decrease rapidly, forming the large peaks in CL and CD 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 CL 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°, CL was approximately 0.6. In fact, from our steady-state calculations and from the measurements of Vogel (1967), the maximum value of CL at steady-state conditions is approximately 0.6. For reference, we used this maximum CL value and the translational velocity of the wing, ut (velocity at r0), to estimate the quasi-steady lift of the wing, 0.6(0.5ρut2S), and therefore the quasi-steady lift coefficient, 0.6(ut2/U2). The quasi-steady lift coefficient is plotted in Fig. 3B for comparison. For the flapping motion described above, CL is much larger than the corresponding quasi-steady value for most of the stroke. The mean value of CL is 1.2, more than twice the quasi-steady value (0.5 in Fig. 3). As seen in Fig. 3, this large mean CL results from three causes: the CL peak at the beginning of the stroke, the large CL in the middle of the stroke and the large CL peak near the end of the stroke. The above analysis showed that the CL peak at the beginning of the stroke could be explained by the rapid acceleration of the wing, the large CL in the middle of the stroke by the stall-absence mechanism and the large CL 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 CL 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 CL peak at the beginning of the stroke, indicating that the CL 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 CL 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 CL, slightly decreasing it instead. This detrimental effect of the wake from the previous stroke on CL 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 CL 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.

### Effects of varying the location of the wing-rotation axis and the rotation rate

The large CL 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 CD peak that accompanies the large CL peak (see Figs 3 and 6) (see also fig. 1D in Dickinson et al., 1999). It was shown above that the large CL and CD 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 CL (and CD) 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 CL peak. These results further show that the large CL peak near the end of a stroke is due to the pitching-up rotation of the wing.

### Effects of the timing of wing rotation

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 CL and CD 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, CL is also smaller (no CL peak exists) because, at this point, the wing is conducting pitching-down rotation, cancelling some of the effect of the translational acceleration. The mean CL 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 CL 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 CL peak near the end of the stroke.

In the rotation-delayed case (Fig. 11), near the end of a stroke, there is no CL peak; instead, CL 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, CL. The mean CL 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 CL 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 CL pattern in the last part of the stroke.

### Comparisons between the model-wing experiment and fruit fly data

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.

#### Comparison with model-wing experiment

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×103 kg m–3, the wing length l was 0.25 m and the wing area S was 0.0167 m2. 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 r0=0.58l) is calculated as U=0.14 m s–1. From the above data,

$$0.5{\rho}\mathit{U}^{2}\mathit{S}{=}0.14{\ }N$$
⁠. Using the definition of CL (equation 23), the lift in fig. 3A of Dickinson et al. (1999) can be converted to CL (Fig. 14). The CL curve from Fig. 6B with Δτt=1.14 (approximately the same as that used by Dickinson et al., 1999) is plotted for comparison. The calculated CL is clearly smaller than the measured values, and they differ by almost the same amount throughout the majority of the stroke. The model wing used in the present calculations had an aspect ratio of 2.76, which is much smaller than that used by Dickinson et al. (1999), which can be calculated as
$$\mathit{l^{2}/S}{=}3.74$$
. This partly explain the lower CL. Despite the quantitative differences between the calculated and experimental results, the overall pattern is very similar.

#### Comparison with fruit-fly 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 St was 0.058 cm2, 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 CL in equation 23:

${\bar{L}}\ {=}\ \overline{C_{L}}\ {\times}\ 0.5{\rho}\mathit{U}^{2}\mathit{S}_{t}{\ },$
25

where is the mean lift and is the mean lift coefficient (air density ρ is taken as 1.25×10–6 g cm–4 s2). The speed at the radial location r0 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φnr0=219 cm s–1, and U can be calculated as U=219/0.825=265 cm s–1. Inserting the values of ρ, U and St into equation 25, the mean lift can be written as

$$\overline{\mathit{L}}{=}\overline{\mathit{C}_{L}}{\times}2.50{\times}10^{{\mbox{--}}5}{\ }N$$
⁠. For to equal the insect weight, should therefore be 1.96×10–5/(2.50×10–5)=0.78.

For the symmetrical rotation and rotation-advanced cases, 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°, was 0.79 for symmetrical rotation and 1.02 for the rotation-advanced case; when the angle of attack was 34°, 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

CD    drag coefficient

CL    lift coefficient

D    drag

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

l    wing length

L    lift

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

r0    radius of the second moment of wing area

Re    Reynolds number

S    area of one wing

St    area of a wing pair

t*    time

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

ut    translational velocity of the wing

ut+    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)

τ1    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

Fig. 1.

Sketches of the reference frames and wing motion. (A) oxyz is an inertial frame, with the xz plane in the stroke plane. oxyz′ is a frame fixed on the wing, with the x′ axis along the wing chord and the z′ axis along the wing span. ψ is the azimuthal angle of the wing, α is the angle of attack and l is the distance between the y axis and the wing tip or the wing length. (B) The motion of a section of the wing.

Fig. 1.

Sketches of the reference frames and wing motion. (A) oxyz is an inertial frame, with the xz plane in the stroke plane. oxyz′ is a frame fixed on the wing, with the x′ axis along the wing chord and the z′ axis along the wing span. ψ is the azimuthal angle of the wing, α is the angle of attack and l is the distance between the y axis and the wing tip or the wing length. (B) The motion of a section of the wing.

Fig. 2.

Portions of the body-conforming grid near the wing surface. (A) In a sectional plane; (B) in the y′=0 plane (see Fig. 1A for a definition of this plane).

Fig. 2.

Portions of the body-conforming grid near the wing surface. (A) In a sectional plane; (B) in the y′=0 plane (see Fig. 1A for a definition of this plane).

Fig. 3.

Wing translational (ut+) and rotational (α̇+) non-dimensional velocities (A), lift coefficient CL (B) and drag coefficient CD (C) plotted against non-dimensional time τ. In B, the quasi-steady lift coefficient is also plotted for comparison.

Fig. 3.

Wing translational (ut+) and rotational (α̇+) non-dimensional velocities (A), lift coefficient CL (B) and drag coefficient CD (C) plotted against non-dimensional time τ. In B, the quasi-steady lift coefficient is also plotted for comparison.

Fig. 4.

(A–H) Vorticity plots at three spanwise locations at various times during one stroke. r1, r2 and r3 denote locations 0.75, 0.5 and 0.25 wing lengths from the wing root, respectively. τ, non-dimensional time; α, angle of attack. Solid and broken vorticity lines indicate positive and negative vorticity, respectively. The magnitude of the non-dimensional vorticity at an outer contour is 1. Starting from the outer contour, for the first 21 contours, the contour interval is 0.2, for the next 30 contours, the contour interval is 0.5, and for the remainder of the contours, the contour interval is 5.

Fig. 4.

(A–H) Vorticity plots at three spanwise locations at various times during one stroke. r1, r2 and r3 denote locations 0.75, 0.5 and 0.25 wing lengths from the wing root, respectively. τ, non-dimensional time; α, angle of attack. Solid and broken vorticity lines indicate positive and negative vorticity, respectively. The magnitude of the non-dimensional vorticity at an outer contour is 1. Starting from the outer contour, for the first 21 contours, the contour interval is 0.2, for the next 30 contours, the contour interval is 0.5, and for the remainder of the contours, the contour interval is 5.

Fig. 5.

Sectional streamline plots at three spanwise locations at various times during one stroke. r1, r2 and r3 denote locations 0.75, 0.5 and 0.25 wing lengths from the wing root, respectively. τ, non-dimensional time; α, angle of attack (the spatial interval of the incoming streamlines can be seen from the right of each plot).

Fig. 5.

Sectional streamline plots at three spanwise locations at various times during one stroke. r1, r2 and r3 denote locations 0.75, 0.5 and 0.25 wing lengths from the wing root, respectively. τ, non-dimensional time; α, angle of attack (the spatial interval of the incoming streamlines can be seen from the right of each plot).

Fig. 6.

Effects of the period of acceleration (Δτt). The wing translational (ut+) and rotational (α̇+) non-dimensional velocities (A), lift coefficient CL (B) and drag coefficient CD (C) plotted against non-dimensional time τ.

Fig. 6.

Effects of the period of acceleration (Δτt). The wing translational (ut+) and rotational (α̇+) non-dimensional velocities (A), lift coefficient CL (B) and drag coefficient CD (C) plotted against non-dimensional time τ.

Fig. 7.

Effects of the wake from previous strokes. Lift CL and drag CD coefficients plotted against non-dimensional time τ during one stroke. The solid lines represent results replotted from Fig. 6; the dashed lines represent results from a wing started in still air. Δτt, period of acceleration.

Fig. 7.

Effects of the wake from previous strokes. Lift CL and drag CD coefficients plotted against non-dimensional time τ during one stroke. The solid lines represent results replotted from Fig. 6; the dashed lines represent results from a wing started in still air. Δτt, period of acceleration.

Fig. 8.

Velocity vectors in the middle section of the wing near the start of a stroke, corresponding to the vorticity field in Fig. 4A. The horizontal arrow at the top right represents the velocity of the wing in the phase of constant-speed translation and it serves as a scale for the velocity vectors in the figure.

Fig. 8.

Velocity vectors in the middle section of the wing near the start of a stroke, corresponding to the vorticity field in Fig. 4A. The horizontal arrow at the top right represents the velocity of the wing in the phase of constant-speed translation and it serves as a scale for the velocity vectors in the figure.

Fig. 9.

Effects of pitching-axis location (as a percentage of mean chord length c from the leading edge of the wing). Lift CL and drag CD coefficients plotted against non-dimensional time τ during one cycle.

Fig. 9.

Effects of pitching-axis location (as a percentage of mean chord length c from the leading edge of the wing). Lift CL and drag CD coefficients plotted against non-dimensional time τ during one cycle.

Fig. 10.

Effects of the pitching-up rotation rate Δτr. The wing translational (ut+) and rotational (α̇+) non-dimensional velocities (A), lift coefficient CL (B) and drag coefficient CD (C) plotted against non-dimensional time τ.

Fig. 10.

Effects of the pitching-up rotation rate Δτr. The wing translational (ut+) and rotational (α̇+) non-dimensional velocities (A), lift coefficient CL (B) and drag coefficient CD (C) plotted against non-dimensional time τ.

Fig. 11.

Effects of the timing of wing rotation. The wing translational (ut+) and rotational (α̇+) non-dimensional velocities (A), lift coefficient CL (B) and drag coefficient CD (C) plotted against non-dimensional time τ. Rotation is defined as advanced, symmetrical or delayed with respect to stroke reversal.

Fig. 11.

Effects of the timing of wing rotation. The wing translational (ut+) and rotational (α̇+) non-dimensional velocities (A), lift coefficient CL (B) and drag coefficient CD (C) plotted against non-dimensional time τ. Rotation is defined as advanced, symmetrical or delayed with respect to stroke reversal.

Fig. 12.

(A–E) Vorticity plots at three spanwise locations at various times during one stroke (symmetrical rotation). r1, r2 and r3 denote locations at 0.75, 0.5 and 0.25 wing lengths from the wing root, respectively. τ, non-dimensional time; α, angle of attack. The magnitude of the non-dimensional vorticity at an outer contour is 1. Starting from the outer contour, for the first 21 contours, the contour interval is 0.2, for the next 30 contours, the contour interval is 0.5, and for the reminder of the contours, the contour interval is 5.

Fig. 12.

(A–E) Vorticity plots at three spanwise locations at various times during one stroke (symmetrical rotation). r1, r2 and r3 denote locations at 0.75, 0.5 and 0.25 wing lengths from the wing root, respectively. τ, non-dimensional time; α, angle of attack. The magnitude of the non-dimensional vorticity at an outer contour is 1. Starting from the outer contour, for the first 21 contours, the contour interval is 0.2, for the next 30 contours, the contour interval is 0.5, and for the reminder of the contours, the contour interval is 5.

Fig. 13.

(A–F) Vorticity plots at three spanwise locations at various times during one stroke (rotation delayed). r1, r2 and r3 denote locations at 0.75, 0.5 and 0.25 wing lengths from the wing root, respectively. τ, non-dimensional time; α, angle of attack. The magnitude of the non-dimensional vorticity at an outer contour is 1. Starting from the outer contour, for the first 21 contours, the contour interval is 0.2, for the next 30 contours, the contour interval is 0.5, and for the reminder of contours, the contour interval is 5.

Fig. 13.

(A–F) Vorticity plots at three spanwise locations at various times during one stroke (rotation delayed). r1, r2 and r3 denote locations at 0.75, 0.5 and 0.25 wing lengths from the wing root, respectively. τ, non-dimensional time; α, angle of attack. The magnitude of the non-dimensional vorticity at an outer contour is 1. Starting from the outer contour, for the first 21 contours, the contour interval is 0.2, for the next 30 contours, the contour interval is 0.5, and for the reminder of contours, the contour interval is 5.

Fig. 14.

Comparison of the lift coefficient CL, calculated here for the rotation-advanced case, with that taken from the experimental data of fig. 3A of Dickinson et al. (1999).

Fig. 14.

Comparison of the lift coefficient CL, calculated here for the rotation-advanced case, with that taken from the experimental data of fig. 3A of Dickinson et al. (1999).

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.

Dickinson, M. H. (
1994
). The effects of wing rotation on unsteady aerodynamic performance at low Reynolds numbers.
J. Exp. Biol
.
192
,
179
–206.
Dickinson, M. H. and Götz, K. G. (
1993
). Unsteady aerodynamic performance of model wings at low Reynolds numbers.
J. Exp. Biol
.
174
,
45
–64.
Dickinson, M. H., Lehman, F. O. and Sane, S. P. (
1999
). Wing rotation and the aerodynamic basis of insect flight.
Science
284
,
1954
–1960.
Ellington, C. P. (
1984
a). The aerodynamics of hovering insect flight. I. The quasi-steady analysis.
Phil. Trans. R. Soc. Lond. B
305
,
1
–15.
Ellington, C. P. (
1984
b). The aerodynamics of hovering insect flight. III. Kinematics.
Phil. Trans. R. Soc. Lond. B
305
,
41
–78.
Ellington, C. P., van den Berg, C. and Willmott, A. P. (
1996
). Leading edge vortices in insect flight.
Nature
384
,
626
–630.
Hilgenstock, A. (
1988
). A fast method for the elliptic generation of three dimensional grid with full boundary control. In Numerical Grid Generation in CFM’88 (ed. S. Sengupta, J. Hauser, P. R. Eiseman and J. F. Thompson), pp. 137–146. Swansea UK: Pineridge Press Ltd.
Lan, S. L. and Sun, M. (
2001
). Aerodynamic properties of a wing performing unsteady rotational motions at low Reynolds number.
Acta Mech
.
149
,
135
–147.
Liu, H., Ellington, C. P., Kawachi, K., Van Den Berg, C. and Willmott, A. P. (
1998
). A computational fluid dynamic study of hawkmoth hovering.
J. Exp. Biol
.
201
,
461
–477.
Rogers, S. E. and Kwak, D. (
1990
). Upwind differencing scheme for the time-accurate incompressible Navier–Stokes equations.
AIAA J
.
28
,
253
–262.
Rogers, S. E., Kwak, D. and Kiris, C. (
1991
). Numerical solution of the incompressible Navier–Stokes equations for steady-state and dependent problems.
AIAA J
.
29
,
603
–610.
Spedding, G. R. (
1992
). The aerodynamics of flight. In Advances in Comparative and Enviromental Physiology, vol. II, Mechanics of Animal Locomotion (ed. R. McN. Alexander), pp. 51–111. London: Springer-Verlag.
Sun, M. and Hamdani, H. (
2001
). High-lift generation by an airfoil performing unsteady motion at low Reynolds number.
Acta Mech. Sinica
17
,
97
–144.
Vogel, S. (
1967
). Flight in Drosophila. III. Aerodynamic characteristics of fly-wing and wing models.
J. Exp. Biol
.
46
,
431
–443.
Weis-Fogh, T. (
1973
). Quick estimates of flight fitness in hovering animals, including novel mechanism for lift production.
J. Exp. Biol
.
59
,
169
–230.