SUMMARY

The unsteady aerodynamic forces of a model fruit fly wing in flapping motion were investigated by numerically solving the Navier–Stokes equations. The flapping motion consisted of translation and rotation [the translation velocity (ut) varied according to the simple harmonic function (SHF), and the rotation was confined to a short period around stroke reversal]. First, it was shown that for a wing of given geometry with ut varying as the SHF, the aerodynamic force coefficients depended only on five non-dimensional parameters, i.e. Reynolds number (Re), stroke amplitude (Φ), mid-stroke angle of attack(αm), non-dimensional duration of wing rotation(Δτr) and rotation timing [the mean translation velocity at radius of the second moment of wing area (U), the mean chord length (c) and c/U were used as reference velocity, length and time, respectively]. Next, the force coefficients were investigated for a case in which typical values of these parameters were used (Re=200;Φ=150°; αm=40°; Δτr was 20%of wingbeat period; rotation was symmetrical). Finally, the effects of varying these parameters on the force coefficients were investigated.

In the Re range considered (20–1800), when Re was above ∼100, the lift(L) and drag(D) coefficients were large and varied only slightly with Re (in agreement with results previously published for revolving wings); the large force coefficients were mainly due to the delayed stall mechanism. However, when Re was below∼100, L decreased and D increased greatly. At such low Re, similar to the case of higher Re, the leading edge vortex existed and attached to the wing in the translatory phase of a half-stroke; but it was very weak and its vorticity rather diffused, resulting in the small L and large D. Comparison of the calculated results with available hovering flight data in eight species(Re ranging from 13 to 1500) showed that when Re was above∼100, lift equal to insect weight could be produced but when Rewas lower than ∼100, additional high-lift mechanisms were needed.

In the range of Re above ∼100, Φ from 90° to 180°and Δτr from 17% to 32% of the stroke period (symmetrical rotation), the force coefficients varied only slightly with Re, Φand Δτr. This meant that the forces were approximately proportional to the square of Φn (n is the wingbeat frequency); thus, changing Φ and/or n could effectively control the magnitude of the total aerodynamic force.

The time course of L(or D) in a half-stroke for ut varying according to the SHF resembled a half sine-wave. It was considerably different from that published previously for ut, varying according to a trapezoidal function (TF) with large accelerations at stroke reversal, which was characterized by large peaks at the beginning and near the end of the half-stroke. However, the mean force coefficients and the mechanical power were not so different between these two cases (e.g. the mean force coefficients for ut varying as the TF were approximately 10% smaller than those for utvarying as the SHF except when wing rotation is delayed).

Introduction

It has been shown that conventional aerodynamic theory, which was based on steady 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, much progress has been made in revealing the unsteady high-lift mechanisms of flapping insect wings.

Dickinson and Götz(1993) measured the aerodynamic forces of an airfoil started rapidly at high angles of attack in the Reynolds number (Re) range of the fruit fly wing (Re=75–225;for a flapping wing, Re is based on the mean chord length and the mean translation velocity at radius of the second moment of wing area). They showed that lift was enhanced by the presence of a dynamic stall vortex, or leading edge vortex (LEV). After the initial start, lift coefficient(CL) of approximately 2 was maintained within 2–3 chord lengths of travel. Afterwards, CL started to decrease due to the shedding of the LEV. But the decrease was not rapid,possibly because the shedding of the LEV was slow at such low Re; and from 3 to 5 chord lengths of travel, CL was still as high as approximately 1.7. The authors considered that because the fly wing typically moved only 2–4 chord lengths each half-stroke, the stall-delaying behavior was more appropriate for models of insect flight than were the steady-state approximations.

Ellington et al. (1996) and van den Berg and Ellington(1997a,b)performed flow-visualization studies on the large hawkmoth Manduca sexta, in tethered forward flight (speed range, 0.4–5.7 m s–1), and on a mechanical model of the hawkmoth wings(Re≈3500). They found that the LEV on the wings did not shed in the translational phases of the half-strokes and that there was a spanwise flow directed from the wing base to the wing tip. Analysis of the momentum imparted to fluid by the vortex wake showed that the LEV could produce enough lift for the weight support. For the hovering case, the hawkmoth wing traveled approximately three chord lengths each half-stroke, whereas for the case of forward flight at high speeds the wing traveled twice as far. The authors suggested that the spanwise flow had prevented the LEV from detaching.

The above studies identified delayed stall as the high-lift mechanism of some small and large insects. Recently, Dickinson et al.(1999) measured the aerodynamic forces on a revolving model fruit fly wing (Re≈75) and showed that stall did not occur and large lift and drag were maintained. Usherwood and Ellington(2002a,b)measured the aerodynamic forces on revolving real and model wings of various insects and a bird (quail) and, for some cases, flow visualization was also conducted. They found that large aerodynamic forces were maintained by the attachment of the LEV for Re≈600 (mayfly) to 15 000 (quail) and for different wing planforms. These results further showed that the delayed stall mechanism was valid for most insects [wing length (R) 2 mm(fruit fly) to 50 mm (hawkmoth)]. The delayed stall mechanism was confirmed by computational fluid dynamics (CFD) analyses(Liu et al., 1998; Wang, 2000; Lan and Sun, 2001).

Dickinson et al. (1999) and Sane and Dickinson (2001), by measuring the aerodynamic forces on a mechanical model of fruit fly wing in flapping motion, showed that when the translation velocity varied according to a trapezoidal function (TF) with large accelerations at stroke reversal and the wing rotation was advanced, in addition to the large forces during the translational phase of a half-stroke, very large force peaks occurred at the beginning and near the end of the half-stroke. Sun and Tang(2002a) and Ramamurti and Sandberg (2002) simulated the flows of model fruit fly wings using the CFD method, based on wing kinematics nearly identical to those used in the experiment of Dickinson et al.(1999). They obtained results qualitatively similar to those of the experiment. The large forces during the translational phase were explained by the delayed stall mechanism(Dickinson et al., 1999). It was suggested that the large force peaks at the beginning of the half-stroke were due to the rapid translational acceleration of the wing and the interaction between the wing and the wake left by the previous strokes(Dickinson et al., 1999; Sun and Tang, 2002a; Birch and Dickinson, 2003), and those near the end of the stroke were due to the effects of wing rotation(Dickinson et al., 1999; Sun and Tang, 2002a).

The experiments on revolving wings (Usherwood and Ellington, 2002a,b; Dickinson et al., 1999) have showed that similar high force coefficients (due to the delayed stall mechanism) are obtained in the Re range of approximately 140 (model fruit fly wing) to 15 000 (quail wing). It is of interest to investigate the aerodynamic force behavior and the delayed stall mechanism at lower Reynolds numbers, because Re for some very small insects is as low as 20(Weis-Fogh, 1973; Ellington, 1984a). In the experiments (Dickinson et al.,1999; Sane and Dickinson,2001) and CFD simulations (Sun and Tang, 2002a; Ramamurti and Sandberg, 2002) on the flapping model fruit fly wings, the translation velocity (ut) of the wing varied as a TF with rapid accelerations at stroke reversal (the stroke positional angle followed a smoothed triangular wave), which was an idealization on the basis of the kinematic data of tethered fruit flies(Dickinson et al., 1999). However, data of free flight in many insects(Ellington, 1984c; Ennos, 1989) showed that ut was close to the simple harmonic function (SHF). Recent data of free-flying fruit fly (Fry et al.,2003) also showed that ut was close to the SHF. When ut varies as the SHF, the time courses of the forces and their sensitivity to some of the kinematic parameters, such as wing-rotation rate and stroke timing, might be significantly different from those when ut varies as a TF with large acceleration at stroke reversal. Therefore, it is of interest to study the aerodynamic forces for the case of ut varying as the SHF.

In the present study, we use the CFD method to simulate the flows of a model insect wing in flapping motion. The flapping motion consists of translation and rotation (Fig. 1). The translation follows the SHF; the rotation is confined to stroke reversal. As will be shown below, for a given wing with ut varying the SHF (in the absence of wing deformation),its aerodynamic force coefficients depend only on the following five non-dimensional parameters: Re, stroke amplitude (Φ), mid-stroke angle of attack (αm), wing-rotation duration(Δτr) and rotation timing (τr). We consider the Re range of 20 to 1800; in addition, we investigate the effects of varying Φ, αm, Δτr andτ r.

Fig. 1.

Sketches of the reference frames and wing motion. OXYZ is an inertial frame, with the XY plane in the horizontal plane. oxyz is a frame fixed on the wing, with the x-axis along the wing chord and the y-axis along the wing span. ϕ, positional angle of the wing; α, geometrical angle of attack of the wing; R, wing length.

Fig. 1.

Sketches of the reference frames and wing motion. OXYZ is an inertial frame, with the XY plane in the horizontal plane. oxyz is a frame fixed on the wing, with the x-axis along the wing chord and the y-axis along the wing span. ϕ, positional angle of the wing; α, geometrical angle of attack of the wing; R, wing length.

Materials and methods

Most of the procedures used in these CFD simulations have been described elsewhere (Sun and Tang, 2002a,b). The planform of the model wing used (Fig. 2) is the same as that of the robotic fruit fly wing used by Dickinson et al. (1999). The wing section is a flat plate of 3% thickness with round leading and trailing edges. The ratio of the wing length (R) to the mean chord length(c) is 3. The radius of the second moment of wing area(r2) is 0.6R (the mean translational velocity at r2 is used as reference velocity in this study). Two coordinate systems are used. One is the inertial coordinate system, OXYZ, and the other is the body-fixed coordinate system, oxyz (Fig. 1).

Fig. 2.

The wing planform used.

Fig. 2.

The wing planform used.

Stroke kinematics

The velocity at the span location r2 due to wing translation is called the translational velocity (ut). ut is assumed to vary as the SHF:
\[\ u_{\mathrm{t}}^{+}=0.5{\pi}\mathrm{sin}(2{\pi}{\tau}{/}{\tau}_{\mathrm{c}}),\]
(1)
where the non-dimensional translational velocity ut+=ut/U(U is the reference velocity); non-dimensional timeτ=tU/c (t is the time); and τcis the non-dimensional period of a wingbeat cycle. The azimuth-rotational speed of the wing is related to ut. Denoting the azimuth-rotational speed as
\({\dot{{\phi}}}\)
, we have
\({\dot{{\phi}}}(t)=u_{t}{/}r_{2}\)
. The geometric angle of attack of the wing is denoted by α. It assumes a constant value except at the start or near the end of a stroke. The constant value is denoted by αm, the mid-stroke angle of attack. Around the stroke reversal, α changes with time and the angular velocity,
\({\dot{{\alpha}}}\)
, is given by:
\[\ {\dot{{\alpha}}}^{+}={\varpi}\{1-\mathrm{cos}[2{\pi}({\tau}-{\tau}_{\mathrm{r}}){/}{\Delta}{\tau}_{\mathrm{r}}]\};{\tau}_{\mathrm{r}}{\leq}{\tau}{\leq}({\tau}_{\mathrm{r}}+{\Delta}{\tau}_{\mathrm{r}}),\]
(2)
where the non-dimensional form
\({\dot{{\alpha}}}^{+}={\dot{{\alpha}}}c{/}U\)
;ϖ is the mean non-dimensional angular velocity of rotation [note thatϖ here is different by a factor of R/r2 from that defined in Ellington (1984c),where velocity at wing tip was used as reference velocity]; τris the non-dimensional time at which the rotation starts;Δτ r is the non-dimensional time interval over which the rotation lasts, which is termed as wing-rotation duration. In the time interval of Δτr, the wing rotates fromα=αm to α=180–αm. Therefore, when αm and Δτr are specified, ϖ can be determined.

In the flapping motion described above, the period of wingbeat cycleτ c, the geometric angle of attack at mid-strokeα m, the rotation duration Δτr or the mean angular velocity rotation ϖ and the rotation timing τrneed to be specified. Note that since U=2Φnr2(where n is the wingbeat frequency and Φ is the stroke amplitude), τc (=U/cn) is related to Φ byτ c=2Φ(r2/R)(R/c).

Flow equations and evaluation of the aerodynamic forces

The governing equations of the flow are the three-dimensional incompressible unsteady Navier–Stokes equations. Written in the inertial coordinate system OXYZ and non-dimensionalized, they are as follows:
\[\ \frac{{\partial}u}{{\partial}X}+\frac{{\partial}{\nu}}{{\partial}Y}+\frac{{\partial}w}{{\partial}Z}=0,\]
(3)
\[\ \frac{{\partial}u}{{\partial}{\tau}}+u\frac{{\partial}u}{{\partial}X}+{\nu}\frac{{\partial}u}{{\partial}Y}+w\frac{{\partial}u}{{\partial}Z}=-\frac{{\partial}p}{{\partial}X}+\frac{1}{Re}\left(\frac{{\partial}^{2}u}{{\partial}X^{2}}+\frac{{\partial}^{2}u}{{\partial}Y^{2}}+\frac{{\partial}^{2}u}{{\partial}Z^{2}}\right),\]
(4)
\[\ \frac{{\partial}{\nu}}{{\partial}{\tau}}+u\frac{{\partial}{\nu}}{{\partial}X}+{\nu}\frac{{\partial}{\nu}}{{\partial}Y}+w\frac{{\partial}{\nu}}{{\partial}Z}=-\frac{{\partial}p}{{\partial}Y}+\frac{1}{Re}\left(\frac{{\partial}^{2}{\nu}}{{\partial}X^{2}}+\frac{{\partial}^{2}{\nu}}{{\partial}Y^{2}}+\frac{{\partial}^{2}{\nu}}{{\partial}Z^{2}}\right),\]
(5)
\[\ \frac{{\partial}w}{{\partial}{\tau}}+u\frac{{\partial}w}{{\partial}X}+{\nu}\frac{{\partial}w}{{\partial}Y}+w\frac{{\partial}w}{{\partial}Z}=-\frac{{\partial}p}{{\partial}Z}+\frac{1}{Re}\left(\frac{{\partial}^{2}w}{{\partial}X^{2}}+\frac{{\partial}^{2}w}{{\partial}Y^{2}}+\frac{{\partial}^{2}w}{{\partial}Z^{2}}\right),\]
(6)
where u, v and w are three components of the non-dimensional fluid velocity and p is the non-dimensional fluid pressure. In the non-dimensionalization, U, c and c/U are taken as reference velocity, length and time, respectively. Re in equations 4–6 is defined as Re=cU/ν (where ν is the kinematic viscosity of the fluid). The numerical method used to solve equations 3–6 is the same as that in Sun and Tang(2002a,b).
Once the Navier–Stokes equations are numerically solved, the fluid velocity components and pressure at discretized grid points for each time step are available. The aerodynamic forces (lift, L, and drag, D)acting on the wing are calculated from the pressure and the viscous stress on the wing surface. The lift and drag coefficients are defined as follows:
\[\ C_{\mathrm{L}}=\frac{L}{0.5{\rho}U^{2}S},\]
(7)
\[\ C_{\mathrm{D}}=\frac{D}{0.5{\rho}U^{2}S},\]
(8)
where ρ is the fluid density and S is the wing area.

Non-dimensional parameters that affect the aerodynamic force coefficients

For a wing of given geometry (in the absence of deformation), when its flapping motion is prescribed, solution of the non-dimensional Navier–Stokes equations (equations 4–6) gives the aerodynamic force coefficients CL and CD; the only non-dimensional parameter in the Navier–Stokes equations that needs to be specified is Re. To prescribe the flapping motion, as mentioned above, Φ, αm, Δτr andτ r need to be specified. That is, the aerodynamic force coefficients on the wing depend on five non-dimensional parameters: Re, Φ, αm, Δτr andτ r. When the wing rotation is symmetrical, τrmay be determined from Δτr; thus, CLand CD depend only on four parameters: Re, Φ,α m and Δτr.

Results

Code validation and grid resolution test

Code validation

The code used in this study is the same as that in Sun and Tang(2002a). It was tested by measured unsteady aerodynamic forces on a flapping model fruit fly wing(Sun and Tang, 2002b; Sun and Wu, 2003). The calculated drag coefficient agreed well with the measured value [see fig. 2A,C of Sun and Wu (2003)]. For the lift coefficient, in the translation phase during the middle, and in the rotation phase at the end, of each half-stroke, the computed value agreed well with the measured value, whereas in the beginning of the stroke, the computed peak value was much smaller than the measured value [see fig. 2B,D of Sun and Wu (2003) and fig. 4 of Sun and Tang (2002b)]. Recently,Birch and Dickinson (2003)visualized the vorticity patterns around the flapping model fruit fly wing using digital particle image velocimetry. It is of interest to compare the vorticity patterns calculated by Sun and Tang(2002a) using the code with the experimentally visualized ones. For convenience, we define a non-dimensional parameter, , such that =0 at the start of the downstroke and =1 at the end of the subsequent upstroke. At the beginning of the half-stroke, difference in the positions of shed vortices exists between the computation and the experiment[compare fig. 4A of Sun and Tang(2002a) with the panel at =0.02 in fig. 5 of Birch and Dickinson (2003)]; during the translation phase at the middle, and the rotation phase at the end, of the half-stroke, the computed vorticity patterns agree well with the experimentally visualized patterns [compare fig. 4B–E,G,H of Sun and Tang (2002a) with panels at =0.07, 0.12, 0.19, 0.26, 0.38 and 0.48 in fig. 5 of Birch and Dickinson(2003)]. The vorticity comparison is consistent with the force comparison described above: both show that discrepancy exists at the beginning of the half-stroke. The discrepancy might be because the CFD code does not resolve satisfactorily the complex flow at stroke reversal. There is also the possibility that it is due to variations in the precise kinematic patterns, especially at stroke reversal.

Upon the suggestion of a referee of the present paper, we made a further test of the code using the recent experimental data of Usherwood and Ellington(2002a,b)on revolving model wings. In the computation, the wing rotated 120° after the initial start, and Re was set as 1800 (this Re value was similar to that of the bumble bee wing). In order to make comparisons with the experimental data, lift and drag coefficients were averaged between 60°and 120° from the end of the initial start of rotation. The computed and measured CL and CD are shown in Fig. 3 [measured data are taken from fig. 7 of Usherwood and Ellington(2002b)]. In the whole αrange (from –20° to 100°), the computed CLagrees well with the measured values; both have approximately sinusoidal dependence on α. The computed CD also agrees well with the measured values except when α is larger than ∼60°.

Fig. 3.

Comparison of the calculated and measured lift (CL) and drag (CD) coefficients. The experimental data are taken from fig. 7 of Usherwood and Ellington(2002b). α, angle of attack.

Fig. 3.

Comparison of the calculated and measured lift (CL) and drag (CD) coefficients. The experimental data are taken from fig. 7 of Usherwood and Ellington(2002b). α, angle of attack.

The above comparisons show that there still exist some discrepancies between the CFD simulations and the experiments but that, in general, the agreement between the computational and experimental aerodynamic forces is good. We think that the present CFD method can calculate the unsteady aerodynamic forces and flows of the model insect wing with reasonable accuracy.

Grid resolution test

Before proceeding to study the physical aspects of the flow, the effects of the grid density, the time step and computational-domain size on the computed solutions were considered. The sensitivity of the computed flow to spatial and time resolution and to the far-field boundary location was evaluated for the case of Re=1800 (this Reynolds number is the highest among the cases considered in this study). Calculations were performed using three different grid systems. Grid 1 had dimensions 53×48×41 (around the wing section, in the normal direction and in the spanwise direction respectively),grids 2 and 3 had dimensions 77×70×61 and 109×93×78,respectively. The spacings at the wall were 0.003, 0.002 and 0.0015 for grids 1, 2 and 3, respectively. The far-field boundary for these three grids was set at 20c away from the wing surface in the normal direction and 8c away from the wing-tips in the spanwise direction. The grid points were clustered densely toward the wing surface and toward the wake.

Fig. 4 shows the time course of the lift coefficient in one cycle and the contours of the non-dimensional spanwise component of vorticity at mid-span location near the end of a half-stroke (just before the wing starting the pitching-up rotation),calculated using the above three grids and a time-step value of 0.02. It is observed that the first grid refinement produced some change in the vorticity plot; however, after the second grid refinement, the discrepancies are considerably reduced. The differences between the computed lift coefficients using the three grids are small; there is almost no difference between the lift coefficients computed using grids 2 and 3. Computations using grid 3 and two time-step values, Δτ=0.02 and 0.01, were conducted. Discrepancies between the computed aerodynamic forces and vorticity fields using the two time steps were very small. Finally, the sensitivity of the solution to the far-field boundary location was considered by calculating the flow in a large computational domain. In order to isolate the effect of the far-field boundary location, the boundary was made further away from the wing by adding more grid points to the normal direction of grid 3. The calculated results showed that there was no need to put the far-field boundary further than that of grid 3. From the above analysis, it was concluded that grid 3 and a time step value of Δτ=0.02 were appropriate for the present study.

Fig. 4.

Effects of grid density on computed lift coefficient(CL) and vorticity field. (A) The time course of CL in one cycle. (B) Vorticity contour plots at half-wing length near the end of a stroke.

Fig. 4.

Effects of grid density on computed lift coefficient(CL) and vorticity field. (A) The time course of CL in one cycle. (B) Vorticity contour plots at half-wing length near the end of a stroke.

Forces and flows of a typical case

We first considered a case in which typical values of wing kinematic parameters were used (Re=200, Φ=150°,α m=40°, Δτr=1.87 and wing rotation was symmetrical; with the above values of Φ, αm andΔτ r, we had τc=9.37, ϖ andΔτ r=0.2τc).

Fig. 5 shows the time courses of CL and CD in one cycle. CL in the middle portion of a half-stroke is large and dominates over CL at the beginning and near the end of the half-stroke (CD behaves similarly). In previous studies by Dickinson et al. (1999) and Sun and Tang (2002a), in which ut varied as a TF with large accelerations at stroke reversal, large force peaks occurred near the end of the half-stroke. They were caused by the pitching-up rotation of the wing while it was still translating at relatively large velocity. In the present case, peaks in CL and CD also exist(Fig. 5B,C) but they are very small. This is because near the end of the half-stroke the ut has become very low and wing rotation cannot produce a large force at low ut.

Fig. 5.

Time history of aerodynamic force coefficients of the typical case. (A)Non-dimensional angular velocity of pitching rotation(

\({\dot{{\alpha}}}^{+}\)
⁠) and azimuthal rotation (
\({\dot{{\phi}}}^{+}\)
); (B) time courses of lift coefficient (CL) and (C) drag coefficient(CD) in one cycle. Reynolds number (Re)=200,stroke amplitude (Φ)=150°, midstroke angle of attack(αm)=40° and non-dimensional duration of wing rotation(Δτr)=1.87 [mean angular velocity of rotation(ϖ)=0.93]; symmetrical rotation.

Fig. 5.

Time history of aerodynamic force coefficients of the typical case. (A)Non-dimensional angular velocity of pitching rotation(

\({\dot{{\alpha}}}^{+}\)
⁠) and azimuthal rotation (
\({\dot{{\phi}}}^{+}\)
); (B) time courses of lift coefficient (CL) and (C) drag coefficient(CD) in one cycle. Reynolds number (Re)=200,stroke amplitude (Φ)=150°, midstroke angle of attack(αm)=40° and non-dimensional duration of wing rotation(Δτr)=1.87 [mean angular velocity of rotation(ϖ)=0.93]; symmetrical rotation.

The mean lift (L) and drag (D) coefficients are 1.66 and 1.67, respectively, which are much larger than the steady-state values [measured steady-state CL and CD on a fruit fly wing in uniform free-stream in a wind tunnel at the same Re (200) and same αm (40°)are 0.6 and 0.75, respectively (Vogel,1967)]. As seen in Fig. 5, the major part of the mean lift (or drag) comes from the mid-portions of the half-strokes. During these periods, the wing is in pure translational motion (α is constant). From the results in Fig. 5, it is estimated that 88% of the mean lift is contributed by the pure translational motion. As was shown previously (Ellington et al.,1996; Liu et al.,1998; Dickinson et al.,1999; Sun and Tang,2002a), the large CL and CD during the translatory phase of a half-stroke were due to the delayed stall mechanism. That is, the delayed stall mechanism is mainly responsible for the large aerodynamic forces produced. The flow-field data provide further evidence for the above statement. The contours of the non-dimensional spanwise component of vorticity at mid-span location are given in Fig. 6. The LEV does not shed in an entire half-stroke, showing that the large CLand CD in the mid-portion of the half-stroke are due to the delayed stall mechanism.

Fig. 6.

Vorticity plots at half-wing length at various times during one cycle. Solid and broken lines indicate positive and negative vorticity, respectively. The magnitude of the non-dimensional vorticity at the outer contour is 2 and the contour interval is 3. A–D, downstroke; E–H, upstroke. Reynolds number (Re)=200, stroke amplitude (Φ)=150°,midstroke angle of attack (αm)=40° and non-dimensional duration of wing rotation (Δτr)=1.87 [mean angular velocity of rotation (ϖ)=0.93]; symmetrical rotation.

Fig. 6.

Vorticity plots at half-wing length at various times during one cycle. Solid and broken lines indicate positive and negative vorticity, respectively. The magnitude of the non-dimensional vorticity at the outer contour is 2 and the contour interval is 3. A–D, downstroke; E–H, upstroke. Reynolds number (Re)=200, stroke amplitude (Φ)=150°,midstroke angle of attack (αm)=40° and non-dimensional duration of wing rotation (Δτr)=1.87 [mean angular velocity of rotation (ϖ)=0.93]; symmetrical rotation.

The effects of Re

Fig. 7 shows the time courses of CL and CD in one cycle for various Re (Re ranging from 20 to 1800; other conditions being the same as in the typical case). In general, CLincreases and CD decreases as Re increases. However, when Re is higher than ∼100, CL and CD do not vary greatly, whereas when Re is lower than ∼100, CL is much smaller and CD much larger than at higher Re.

Fig. 7.

Time courses of the (A) lift coefficient (CL) and (B)drag coefficient (CD) in one cycle for various Reynolds number (Re). Stroke amplitude (Φ)=150°, midstroke angle of attack (αm)=40° and non-dimensional duration of wing rotation (Δτr)=1.87 [mean angular velocity of rotation(ϖ)=0.93]; symmetrical rotation.

Fig. 7.

Time courses of the (A) lift coefficient (CL) and (B)drag coefficient (CD) in one cycle for various Reynolds number (Re). Stroke amplitude (Φ)=150°, midstroke angle of attack (αm)=40° and non-dimensional duration of wing rotation (Δτr)=1.87 [mean angular velocity of rotation(ϖ)=0.93]; symmetrical rotation.

For the case of Re=200, as discussed above, the large CL and CD during a stroke are due to the delayed stall mechanism. For the cases of other Re, as seen in Fig. 7, CLand CD do not have a sudden drop during a half-stroke(between τ=0 and τ=0.5τc, the downstroke; betweenτ=0.5τc and τ=τc, the upstroke), i.e. stall is also delayed for an entire half-stroke. Fig. 8 shows the vorticity contour plots at half-wing length near the end of a half-stroke. It is seen that, at all Re considered, the LEV does not shed and the delayed stall mechanism exists. However, for Re lower than ∼100, the LEV is very diffused and weak compared with that for higher Re (comparing Fig. 8D,E with Fig. 8A–C; the strength of the LEV can be estimated from the values of vorticity represented by the contours and the spacing between the contours), resulting in small CL and large CD. For reference,vorticity contour plots at various times in one cycle for the case of Re=20 are shown in Fig. 9.

Fig. 8.

Vorticity plots at half-wing length near the end of a half-stroke at various Reynolds number (Re). Solid and broken lines indicate positive and negative vorticity, respectively. The magnitude of the non-dimensional vorticity at the outer contour is 2 and the contour interval is 3. Stroke amplitude (Φ)=150°, midstroke angle of attack(αm)=40° and non-dimensional duration of wing rotation(Δτr)=1.87 [mean angular velocity of rotation(ϖ)=0.93]; symmetrical rotation.

Fig. 8.

Vorticity plots at half-wing length near the end of a half-stroke at various Reynolds number (Re). Solid and broken lines indicate positive and negative vorticity, respectively. The magnitude of the non-dimensional vorticity at the outer contour is 2 and the contour interval is 3. Stroke amplitude (Φ)=150°, midstroke angle of attack(αm)=40° and non-dimensional duration of wing rotation(Δτr)=1.87 [mean angular velocity of rotation(ϖ)=0.93]; symmetrical rotation.

Fig. 9.

Vorticity plots at half-wing length at various times during one cycle. Solid and broken lines indicate positive and negative vorticity, respectively. The magnitude of the non-dimensional vorticity at the outer contour is 2 and the contour interval is 3. A-D, downstroke; E-H, upstroke. Reynolds number Re=20, stroke amplitude Φ=150°, midstroke angle of attackα m=40° and non-dimensional duration of wing rotationΔτ r=1.87 (mean angular velocity of rotation ϖ=0.93);symmetrical rotation.

Fig. 9.

Vorticity plots at half-wing length at various times during one cycle. Solid and broken lines indicate positive and negative vorticity, respectively. The magnitude of the non-dimensional vorticity at the outer contour is 2 and the contour interval is 3. A-D, downstroke; E-H, upstroke. Reynolds number Re=20, stroke amplitude Φ=150°, midstroke angle of attackα m=40° and non-dimensional duration of wing rotationΔτ r=1.87 (mean angular velocity of rotation ϖ=0.93);symmetrical rotation.

L and D at various Reare plotted in Fig. 10. For Re above ∼100, change in L and D with Re are small, whereas for Re below ∼100, L decreases and D increases rapidly as Re decreases. Similar to the typical case, approximately 85–90%of L is contributed by the pure translational motion for all the values of Re considered.

Fig. 10.

Mean lift (L) and drag(D) coefficients vs Reynolds number (Re). Stroke amplitude (Φ)=150°,mid-stroke angle of attack (αm)=40° and non-dimensional duration of wing rotation (Δτr)=1.87 [mean angular velocity of rotation (ϖ)=0.93]; symmetrical rotation.

Fig. 10.

Mean lift (L) and drag(D) coefficients vs Reynolds number (Re). Stroke amplitude (Φ)=150°,mid-stroke angle of attack (αm)=40° and non-dimensional duration of wing rotation (Δτr)=1.87 [mean angular velocity of rotation (ϖ)=0.93]; symmetrical rotation.

Dependence of the force coefficients on mid-stroke angle of attack

Fig. 11 gives L and D in the range ofα m from 25° to 60° (for all Re considered in the above section). The slope of the Lm)curve is approximately constant between αm=25° and 35°; beyond αm=35°, it decreases gradually to zero atα m≈50°.

Fig. 11.

Mean lift (L) and drag(D) coefficients vs mid-stroke angle of attack (αm) for various Reynolds number (Re). Stroke amplitude (Φ)=150° and non-dimensional duration of wing rotation (Δτr)=1.87[mean angular velocity of rotation (ϖ)=0.93]; symmetrical rotation.

Fig. 11.

Mean lift (L) and drag(D) coefficients vs mid-stroke angle of attack (αm) for various Reynolds number (Re). Stroke amplitude (Φ)=150° and non-dimensional duration of wing rotation (Δτr)=1.87[mean angular velocity of rotation (ϖ)=0.93]; symmetrical rotation.

The rate of change of Lwith αm(dL/dαm)from αm=25° to 35° is given in Table 1. For Re above∼100,dL/dαmhardly varies with Re and its value is approximately 3.0, which is almost the same as the measured value (2.9–3.1) for the revolving wings[see fig. 6 of Usherwood and Ellington(2002b); the cited value is for the case of aspect ratio equal to 6 (R/c=3)]. For Rebelow ∼100,dL/dαmdecreases greatly.

Table 1.

The rate of change of lift coefficient(L) withα m(dL/dαm)from αm=25° to 35°

Re 1800 600 200 60 20 
dL/dαm 3.0 2.9 2.9 2.5 1.8 
Re 1800 600 200 60 20 
dL/dαm 3.0 2.9 2.9 2.5 1.8 

Stroke amplitude (Φ)=150°; non-dimensional rotation duration(Δτr)=1.87 (=0.2τc; τc,non-dimensional wingbeat period); symmetrical rotation.

The effects of rotation duration

In the calculations above, Δτr=1.87(=0.2τc; ϖ=0.93). Observation of many insects in free flight (Ellington, 1984c; Ennos, 1989) showed that ϖranged approximately from 0.8 to 1.4. Here, we investigate the effects of varying Δτr (i.e. varying ϖ) on the aerodynamic force coefficients.

Fig. 12 gives the time courses of CL and CD in one cycle for four values of Δτr; Table 2 gives the mean force coefficients. Varying Δτr does not change the mean force coefficients greatly (see Table 2); when Δτr is almost doubled (varied from 1.27 to 2.40), L and D change only approximately 3%. CL and CD in the mid-portion of a half-stroke vary little with Δτr (see Fig. 12). The force peaks around the stroke reversal are due to the effects of wing rotation(Dickinson et al., 1999; Sun and Tang, 2002a); at a given translation velocity, the peaks increase with rotation rate(Sane and Dickinson, 2002; Hamdani and Sun, 2000). WhenΔτ r is relatively short (ϖ is relatively large), the force peaks are relatively large but they occupy a short period; whenΔτ r is longer (ϖ is smaller), the force peaks become smaller but they occupy a longer period. As a result, the force peaks around the stroke reversal for the cases of different Δτr give more or less the same contribution to the corresponding mean force coefficient. This explains why L and D do not change greatly with Δτr (or ϖ).

Fig. 12.

The effects of rotation duration (Δτr) on force coefficients. (A) Non-dimensional angular velocity of pitching rotation(

\({\dot{{\alpha}}}^{+}\)
⁠) and azimuthal rotation (
\({\dot{{\phi}}}^{+}\)
); (B) time courses of lift coefficient (CL) and (C) drag coefficient(CD) in one cycle. Reynolds number (Re)=200,stroke amplitude (Φ)=150° and mid-stroke angle of attack(αm)=40°; symmetrical rotation.

Fig. 12.

The effects of rotation duration (Δτr) on force coefficients. (A) Non-dimensional angular velocity of pitching rotation(

\({\dot{{\alpha}}}^{+}\)
⁠) and azimuthal rotation (
\({\dot{{\phi}}}^{+}\)
); (B) time courses of lift coefficient (CL) and (C) drag coefficient(CD) in one cycle. Reynolds number (Re)=200,stroke amplitude (Φ)=150° and mid-stroke angle of attack(αm)=40°; symmetrical rotation.

Table 2.

Effects of non-dimensional rotation duration (Δτr)on mean lift (L) and drag(D) coefficients

Δτr
\({\bar{{\omega}}}\)
LD
2.39 0.73 1.60 1.71 
1.88 0.93 1.61 1.70 
1.54 1.13 1.63 1.72 
1.31 1.33 1.65 1.76 
Δτr
\({\bar{{\omega}}}\)
LD
2.39 0.73 1.60 1.71 
1.88 0.93 1.61 1.70 
1.54 1.13 1.63 1.72 
1.31 1.33 1.65 1.76 

Reynolds number (Re)=200; mid-stroke angle of attack(αm)=40°; stroke amplitude (Φ)=120°; symmetrical rotation.

\({\bar{{\omega}}}\)
mean non-dimensional angular velocity of rotation.

The effects of rotation timing

Fig. 13 shows the time courses of CL and CD in one cycle for different rotation timing (τr can be read from Fig. 13A; Rem, Φ and Δτr are the same as those in the typical case). In the case of advanced rotation (the major part of rotation is conducted before stroke reversal), the peaks in CL and CD near the end of a half-stroke are larger than those in the case of symmetrical rotation; this is because the wing conducts pitching-up rotation at a higher translational velocity (see Fig. 13A). At the beginning of the next half-stroke, CL and CD are also larger than their counterparts in the case of symmetrical rotation; this is because the wing does not conduct pitching-down rotation in this period (the wing rotation is almost finished before this period). In the case of delayed rotation (the major part of rotation is conducted after stroke reversal), no CL and CD peaks occur near the end of the half-stroke because the wing does not rotate in this period; in the beginning of the next half-stroke, CL is negative and CD is large compared with that in the case of symmetrical rotation because all of the wing rotation is conducted in this period and the rotation is pitching-down rotation.

Fig. 13.

The effects of rotation timing on force coefficients. (A) Non-dimensional angular velocity of pitching rotation(

\({\dot{{\alpha}}}^{+}\)
⁠) and azimuthal rotation (
\({\dot{{\phi}}}^{+}\)
); (B) time courses of lift coefficient (CL) and (C) drag coefficient(CD) in one cycle. Reynolds number (Re)=200,stroke amplitude (Φ)=150°, mid-stroke angle of attack(αm)=40° and non-dimensional duration of wing rotation(Δτr)=1.87 [mean angular velocity of rotation(ϖ)=0.93].

Fig. 13.

The effects of rotation timing on force coefficients. (A) Non-dimensional angular velocity of pitching rotation(

\({\dot{{\alpha}}}^{+}\)
⁠) and azimuthal rotation (
\({\dot{{\phi}}}^{+}\)
); (B) time courses of lift coefficient (CL) and (C) drag coefficient(CD) in one cycle. Reynolds number (Re)=200,stroke amplitude (Φ)=150°, mid-stroke angle of attack(αm)=40° and non-dimensional duration of wing rotation(Δτr)=1.87 [mean angular velocity of rotation(ϖ)=0.93].

The mean force coefficients are given in Table 3. L and D for the case of advanced rotation are approximately 40% and 30% larger than those for the case of delayed rotation, respectively.

Table 3.

Effects of rotation timing on mean lift(L) and drag(D) coefficients

Rotation timingLDL/D
Symmetrical 1.66 1.67 0.99 
Advanced 1.84 2.11 0.87 
Delayed 1.32 1.61 0.82 
Rotation timingLDL/D
Symmetrical 1.66 1.67 0.99 
Advanced 1.84 2.11 0.87 
Delayed 1.32 1.61 0.82 

Reynolds number (Re)=200; mid-stroke angle of attack(αm)=40°, stroke amplitude (Φ)=150°;non-dimensional rotation duration (Δτr)=1.87(=0.2τc; τc, non-dimensional wingbeat period).

The effects of stroke amplitude

Free-flight data collected from many insects(Ellington, 1984c; Ennos, 1989; Fry et al., 2003) showed that the stroke amplitude, Φ, ranged approximately from 90° to 180°. Moreover, an insect might change Φ to control its aerodynamic force (e.g. Ellington, 1984c; Lehmann and Dickinson, 1998). Here, we investigate the effects of Φ on the force coefficients. Calculations were made for various Φ (65°, 90°, 120°, 150°and 180°) while other parameters were fixed (they are the same as those in the typical case). Fig. 14shows the time courses of CL and CD in one cycle; Table 4 gives the mean force coefficients. Note that when Φ is varied, the non-dimensional period of wingbeat cycle will change(τc=2Φr2/c); sinceΔτ r (i.e. ϖ) is fixed,Δτ rc is different for different Φ. In the range of Φ from 90° to 180°, the effects of varying Φon the force coefficients are not large; when Φ increases or decreases by 30°, L and D change less than 3% and 6%, respectively. When Φ is below approximately 90°, the effects of varying Φ become larger (see the results for Φ=65°; Fig. 14; Table 4). It is of interest to point out the fact that Land D hardly vary withΦ (in the range of Φ from 90° to 180°) means that the mean lift and mean drag vary as Φ2, because the forces are non-dimensionalized by U2, and U equals 2Φnr2 (n is the wingbeat frequency).

Fig. 14.

The effects of stroke amplitude (Φ) on force coefficients. (A)Non-dimensional angular velocity of pitching rotation(

\({\dot{{\alpha}}}^{+}\)
⁠) and azimuthal rotation (
\({\dot{{\phi}}}^{+}\)
); (B) time courses of lift coefficient (CL) and (C) drag coefficient(CD) in one cycle. Reynolds number (Re)=200,mid-stroke angle of attack (αm)=40° and non-dimensional duration of wing rotation (Δτr)=1.87 [mean angular velocity of rotation (ϖ)=0.93]; symmetrical rotation.

Fig. 14.

The effects of stroke amplitude (Φ) on force coefficients. (A)Non-dimensional angular velocity of pitching rotation(

\({\dot{{\alpha}}}^{+}\)
⁠) and azimuthal rotation (
\({\dot{{\phi}}}^{+}\)
); (B) time courses of lift coefficient (CL) and (C) drag coefficient(CD) in one cycle. Reynolds number (Re)=200,mid-stroke angle of attack (αm)=40° and non-dimensional duration of wing rotation (Δτr)=1.87 [mean angular velocity of rotation (ϖ)=0.93]; symmetrical rotation.

Table 4.

Effects of stroke amplitude (Φ) on mean lift(L) and drag(D)coefficients

ΦτcLD
180° 11.24 1.69 1.65 
150° 9.37 1.66 1.67 
120° 7.50 1.61 1.70 
90° 5.62 1.57 1.80 
65° 4.06 1.58 2.08 
ΦτcLD
180° 11.24 1.69 1.65 
150° 9.37 1.66 1.67 
120° 7.50 1.61 1.70 
90° 5.62 1.57 1.80 
65° 4.06 1.58 2.08 

Reynolds number (Re)=200, mid-stroke angle of attack(αm)=40°; non-dimensional rotation duration(Δτr)=1.87; symmetrical rotation. τc,non-dimensional wingbeat period.

Sane and Dickinson (2001)studied the effects of varying Φ and other parameters using a dynamically scaled mechanical model of the fruit fly. Their results [see fig. 5A,C of Sane and Dickinson (2001)] showed that when Φ fell below approximately 120°, L decreased and D increased with decreasing Φ (Dincreased rapidly as Φ became small). In the present simulation(Fig. 14; Table 4), when Φ is below∼90°, we also found D increasing and L decreasing with a decrease in Φ, but the rates of change in L and D are smaller than those reported by Sane and Dickinson(2001). In their experiment,when Φ changed, Re and Δτr also changed, but the ratio of Δτrc did not change; in the present simulation, Re and Δτr did not change when Φ changed. To make further comparison with their results, we made some calculations in which Re and Δτr changed with Φ but Δτrc was kept unchanged(=0.2). The results are given in Fig. 15 and Table 5. The trends of variation in Land D with Φ are similar to those in Sane and Dickinson(2001): when Φ is below approximately 120°, Ldecreases and D increases with Φ decreasing, and when Φ is below 90°, D increases rapidly.

Fig. 15.

The effects of stroke amplitude (Φ) on force coefficients whenΔτ rc is fixed (=20%). (A)Non-dimensional angular velocity of pitching rotation(

\({\dot{{\alpha}}}^{+}\)
⁠) and azimuthal rotation (
\({\dot{{\phi}}}^{+}\)
); (B) time courses of lift coefficient (CL) and (C) drag coefficient(CD) in one cycle. Mid-stroke angle of attack(αm)=40°; symmetrical rotation. Δτr,non-dimensional duration of wing rotation; τc, non-dimensional wingbeat period.

Fig. 15.

The effects of stroke amplitude (Φ) on force coefficients whenΔτ rc is fixed (=20%). (A)Non-dimensional angular velocity of pitching rotation(

\({\dot{{\alpha}}}^{+}\)
⁠) and azimuthal rotation (
\({\dot{{\phi}}}^{+}\)
); (B) time courses of lift coefficient (CL) and (C) drag coefficient(CD) in one cycle. Mid-stroke angle of attack(αm)=40°; symmetrical rotation. Δτr,non-dimensional duration of wing rotation; τc, non-dimensional wingbeat period.

Table 5.

Effects of stroke amplitude (Φ) on mean lift(L) and drag(D) coefficients when Δτrc is fixed (=20%)

ΦτcReΔτrLD
180° 11.24 240 2.25 1.71 1.63 
150° 9.37 200 1.87 1.66 1.67 
120° 7.50 160 1.50 1.59 1.76 
90° 5.62 120 1.12 1.49 1.94 
60° 3.75 80 0.75 1.39 2.52 
ΦτcReΔτrLD
180° 11.24 240 2.25 1.71 1.63 
150° 9.37 200 1.87 1.66 1.67 
120° 7.50 160 1.50 1.59 1.76 
90° 5.62 120 1.12 1.49 1.94 
60° 3.75 80 0.75 1.39 2.52 

Mid-stroke angle of attack (αm)=40°; symmetrical rotation. Re, Reynolds number; Φ, stroke amplitude;τ c, non-dimensional wingbeat period; Δτr,non-dimensional rotation duration.

Discussion

The influence of Re and comparison between the lift coefficients and insect flight data

Previous studies on revolving wings (Usherwood and Ellington, 2002a,b; Dickinson et al., 1999) showed that large aerodynamic force coefficients were produced due to the delayed stall mechanism in the Re range of approximately 140 (model fruit fly wing) to 15 000 (quail wing) and that the force coefficients were not sensitive to Re. The present study on a flapping wing has provided results for lower Re. As seen in Fig. 10, when Re is above ∼100, L and D vary only slightly with Re, in agreement with the previous results. However, when Reis below ∼100, Ldecreases and D increases greatly. This is because at such low Re (20, 60), although the LEV still exists and attaches to the wing in the translational phases of the half-strokes, it is rather weak and its vorticity is considerably diffused(see Figs 8D,E, 9).

From the flight data of an insect, the mean lift coefficient needed for supporting its weight (denoted by L,W) can be determined. Data of free hovering (or very low-speed) flight in eight species were obtained. Six species were from Ellington(1984b,c)[the wing length of these species ranges from 9.3 mm (in Episyrphus balteatus) to 14.1 mm (in Bombus hortorum)]; two smaller ones, Drosophila virilis and Encarsiaformosa, were from Weis-Fogh (1973). These data include: insect mass (M), wing length, mean chord length, radius of second moment of wing area, stroke amplitude and wingbeat frequency (see Table 6). On the basis of these data, the reference velocity, Re and mean lift coefficient needed for supporting insect weight were computed(U=2Φnr2, Re=Uc/ν and L,W=mg/0.5ρU2St,where g and St were the gravitational acceleration and the area of both wings, respectively). Re and L,W are given in Table 6.

Table 6.

Data of free hovering flight and lift coefficient needed for weight support

SpeciesM (mg)R (mm)c (mm)r2/RΦ (deg.)n (s-1)ReL,W
Coleoptera: beetles         
    Coccinella 7-punctata 34.4 11.2 3.23 0.53 177 54 443 1.82 
Diptera: flies, mosquitoes         
    Drosophila virilis 2.0 3.0 0.97 0.58 150 240 147 1.15 
    Tipala obsolete 11.4 12.7 2.38 0.6 123 45.5 245 1.31 
    Episyrphus balteatus 27.3 9.3 2.20 0.57 90 160 408 1.52 
    Episyrphus tenax 68.4 11.4 3.19 0.53 109 157 812 1.10 
Hymenoptera: bees and wasps         
    Apis mellifera 101.9 9.8 3.08 0.54 131 197 1018 1.19 
    Bombus hortorum 226 14.1 4.2 0.54 120 152 1463 1.21 
    Encarsia formosa 0.025 0.62 0.23 0.69 135 400 13 2.87 
    Encarsia formosa0.025 0.65 0.38 0.69 135 400 22 1.62 
SpeciesM (mg)R (mm)c (mm)r2/RΦ (deg.)n (s-1)ReL,W
Coleoptera: beetles         
    Coccinella 7-punctata 34.4 11.2 3.23 0.53 177 54 443 1.82 
Diptera: flies, mosquitoes         
    Drosophila virilis 2.0 3.0 0.97 0.58 150 240 147 1.15 
    Tipala obsolete 11.4 12.7 2.38 0.6 123 45.5 245 1.31 
    Episyrphus balteatus 27.3 9.3 2.20 0.57 90 160 408 1.52 
    Episyrphus tenax 68.4 11.4 3.19 0.53 109 157 812 1.10 
Hymenoptera: bees and wasps         
    Apis mellifera 101.9 9.8 3.08 0.54 131 197 1018 1.19 
    Bombus hortorum 226 14.1 4.2 0.54 120 152 1463 1.21 
    Encarsia formosa 0.025 0.62 0.23 0.69 135 400 13 2.87 
    Encarsia formosa0.025 0.65 0.38 0.69 135 400 22 1.62 

M, insect mass; R, wing length; c, mean chord length; r2, radius of second moment of wing area; Φ,stroke amplitude; n, wingbeat frequency; Re, Reynolds number; L,W, mean lift coefficient needed to balance insect weight. Data of Drosophila virilis and Encarsia formosa are from Weis-Fogh(1973); data of other species are from Ellington(1984b,c).*Wing area extended to include the brim hairs(Ellington, 1975).

Now, we compare the data in Table 6 with the results of model-wing simulation in Fig. 11 [here, we assume that the wing planform does not have a significant effect on lift coefficient; this is true for revolving wings (Usherwood and Ellington, 2002b)]. Of the insects considered, Encarsia formosa has the lowest Re (13) and its L,W is 2.87; when its wing area is extended to include the brim hairs(Ellington, 1975), its L,W is still as high as 1.62. As seen in Fig. 11, at such a low Re, the maximum L is ∼1.15 (atα m≈45°), which is much smaller than its L,W. In the computations that gave the results in Fig. 11, symmetrical rotation was used and Φ was 150°. For reference, we made another calculation in which advanced rotation was used,Φ=180° and αm=45° (this combination of parameters was expected to maximize L). The computation gave L=1.25, which was also much smaller than the L,Wof Encarsia formosa. These results show that using the flapping motion described above, the insect could not produce enough lift to support its weight; i.e. at such low Re, high-lift mechanisms, in addition to the delayed stall mechanism, are needed [Weis-Fogh(1973) suggested the `clap and fling' mechanism]. For other insects, Re is above 100 and, as seen in Fig. 11, at anα m between 30° and 50°, a L equal to L,W can be produced.

The above comparison shows that when Re is higher than ∼100,the delayed stall mechanism can produce enough lift for supporting the insect's weight and when Re is lower than ∼100, additional high-lift mechanisms are needed.

Lift and drag vary approximately with the square of Φn

The non-dimensional Navier–Stokes equations (equations 3–6),the equations prescribing the flapping motion (equations 1, 2) and the equations defining the aerodynamic force coefficients (equations 7, 8) show that the mean force coefficients of a wing of given geometry with ut varying as the SHF depend only on Rem, Φ and Δτr (assuming symmetrical rotation). As already discussed above, when Re is above ∼100, the force coefficients vary only slightly with Re; results in Tables 2, 4 show that the force coefficients vary only slightly with Δτr and also vary only slightly with Φ in the range of Φ approximately from 90° to 180°.

When Φ and/or n is varied, Re will change (note that Re=2Φnr2c/ν). Since the force coefficients hardly vary with Φ and Re, the mean lift() and drag() vary approximately with(Φn)2. That is, changing Φ and/or n can effectively control the aerodynamic forces. For instance, increasing Φ by 15%, can be increased by approximately 32% [note that by increasing αm by 15% (e.g. from 40° to 46°), increases only by ∼10% (see Fig. 11)].

In the above discussion, we have assumed that the force coefficients hardly vary with Re, Φ and Δτr (or ϖ). However,in testing the effects of a particular parameter on the force coefficients, we varied one parameter while keeping all others the same as in the typical case. When the parameters are simultaneously varied, do the force coefficients still vary only slightly? We conducted some further calculations in which, for a range of αm from 25° to 60°, Re, Φ andϖ were simultaneously increased by 20% from those of the typical case. Fig. 16 shows the results. At all αm considered, the force coefficients vary only slightly.

Fig. 16.

Mean lift (L) and drag(D) coefficients vs mid-stroke angle of attack (αm; symmetrical rotation). Re, Reynolds number; Φ, stroke amplitude; ϖ, mean non-dimensional angular velocity of wing rotation.

Fig. 16.

Mean lift (L) and drag(D) coefficients vs mid-stroke angle of attack (αm; symmetrical rotation). Re, Reynolds number; Φ, stroke amplitude; ϖ, mean non-dimensional angular velocity of wing rotation.

Comparison of the present results with those of ut varying as the TF

In some recent experimental (Dickinson et al., 1999; Sane and Dickinson, 2001) and computational (Sun and Tang, 2002a,b; Ramamurti and Sandberg, 2002; Sun and Wu, 2003) studies, ut varying as a TF with rapid accelerations at stroke reversal has been employed. It is of interest to discuss the differences between the present results for ut varying as the SHF and those for ut varying as a TF with rapid accelerations at stroke reversal.

The force coefficients

The time courses of force coefficients are considerably different between the two cases; for the case of ut varying as a TF with rapid accelerations at stroke reversal, the CL (or CD) curve is flat in the mid-portion of a half-stroke and has large peaks at the beginning and near the end of the half-stroke [see fig. 3A,B of Dickinson et al. (1999)and fig. 6 of Sun and Tang(2002a)], whereas for the case of ut varying as the SHF, the CL (or CD) curve grossly resembles a half sine-wave (see Fig. 5). As a result, very large time gradients of aerodynamic force exist in each half-stroke in the case of ut varying as the TF but not in the case of ut varying as the SHF.

However, in spite of the large differences in instantaneous force coefficients, the mean force coefficients are not so different. To examine the quantitative differences of the mean force coefficients between the two cases,we made two sets of computations. In the first set, utvaried as a TF with rapid accelerations at stroke reversal [the duration of translational acceleration at stroke reversal was 0.18τc,similar to that used in Dickinson et al.(1999) and Sun and Tang(2002a,b)];symmetrical rotation, advanced rotation (the major part of rotation conducted before stroke reversal) and delayed rotation (the major part of rotation conducted after the stroke reversal) were considered; other conditions(Re, Φ, αm, Δτr) were the same as those in the typical case. In the second set, utwas replaced by the SHF. The results showed that when wing rotation was symmetrical or advanced, Land D for ut varying as the TF are approximately 12% and 8% smaller than those for ut varying as the SHF, respectively, and when wing rotation was delayed, L for ut varying as the TF was approximately 35% smaller than that for ut varying as the SHF but D was approximately the same for the two cases.

Power requirements

In the studies on power requirements of fruit fly flight by Sun and Tang(2002b) and Sun and Wu(2003), a TF similar to that used in Sun and Tang (2002a)was used for ut (wing rotation was symmetrical). As seen above, L and D for ut varying as the TF are 12% and 8% smaller, respectively,than those for ut varying as the SHF. It is of interest to know the effects of replacing the TF by the SHF on the results presented in Sun and Tang (2002a) and Sun and Wu (2003).

To quantify the effects of the new kinematic model, we made calculations in which the same model wing and kinematic parameters as those in Sun and Tang(2002b) and Sun and Wu(2003) were used, except that the TF for ut was replaced by the SHF. For hover flight,the results in Sun and Tang(2002b) and Sun and Wu(2003) are as follows: mean lift equal to the insect weight is produced at αm=36.5°and the body-mass-specific power is 29 W kg–1; with ut varying as the SHF, the correspondingα m and body-mass-specific power are 30.5° and 31.5 W kg–1, respectively. That is, with the SHF, theα m needed is a few degrees smaller and the body-mass-specific power is approximately 10% larger than that with the TF (similar results were obtained for forward flight).

The reason for the needed αm becoming smaller is obvious. Seeing that L and D with the SHF are larger than their counterparts with the TF by approximately the same percentage (i.e. the L to D ratio is not very different for the two cases), one might expect that the power results with the SHF are approximately the same as those with the TF. However, as seen above,the specific power becomes a little larger. This is because in the case of ut varying as the SHF, both CD and ut at the middle portion of a stroke are larger than their counterparts in the case of ut varying as the TF (note that aerodynamic power is proportional to the mean of the product of CD and ut over a stroke cycle, not to D).

List of symbols

     
  • c

    mean chord length

  •  
  • CD

    drag coefficient

  •  
  • D

    mean drag coefficient

  •  
  • CL

    lift coefficient

  •  
  • L

    mean lift coefficient

  •  
  • L,W

    mean lift coefficient for supporting the insect's weight

  •  
  • D

    drag

  •  
  • mean drag

  •  
  • L

    lift

  •  
  • mean lift

  •  
  • M

    mass of insect

  •  
  • n

    wingbeat frequency

     
  • O,o

    origins of the inertial frame of reference and the non-inertial frame of reference

  •  
  • p

    non-dimensional fluid pressure

  •  
  • R

    wing length

  •  
  • r2

    radius of the second moment of wing area

  •  
  • Re

    Reynolds number

  •  
  • S

    area of one wing

  •  
  • St

    area of a wing pair

  •  
  • t

    time

  •  
  • non-dimensional parameter (being zero at the start of downstroke and 1 at the end of the subsequent upstroke)

  •  
  • U

    reference velocity

  •  
  • ut

    translational velocity of the wing

  •  
  • ut+

    non-dimensional translational velocity of the wing

  •  
  • x,y,z

    coordinates in non-inertial frame of reference

  •  
  • X,Y,Z

    coordinates in inertial frame of reference (Z in vertical direction)

  •  
  • α

    geometric angle of attack

  •  
  • \({\dot{{\alpha}}}\)

    angular velocity of pitching rotation

  •  
  • \({\dot{{\alpha}}}^{+}\)

    non-dimensional angular velocity of pitching rotation

  •  
  • αm

    mid-stroke geometric angle of attack

  •  
  • ϕ

    azimuthal or positional angle

  •  
  • \({\dot{{\phi}}}\)

    angular velocity of azimuthal rotation

  •  
  • \({\dot{{\phi}}}^{+}\)

    non-dimensional angular velocity of azimuthal rotation

  •  
  • Φ

    stroke amplitude

  •  
  • ν

    kinematic viscosity

  •  
  • ρ

    density of fluid

  •  
  • τ

    non-dimensional time

  •  
  • τc

    period of one wingbeat cycle (non-dimensional)

  •  
  • τr

    time when pitching rotation starts (non-dimensional)

  •  
  • Δτr

    duration of wing rotation or flip duration (non-dimensional)

  •  
  • ϖ

    mean non-dimensional angular velocity of wing rotation

Acknowledgements

We thank the two referees whose helpful comments and valuable suggestions greatly improved the quality of the paper. This research was supported by the National Natural Science Foundation of China (10232010).

References

Birch, J. M. and Dickinson, M. H. (
2003
). The influence of wing–wake interactions on the production of aerodynamic forces in flapping flight.
J. Exp. Biol.
206
,
2257
-2272.
Dickinson, M. H. and Gö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. (
1975
). Non-steady-state aerodynamics of the flight of Encarsia formosa. In
Swimming and Flying in Nature
, vol.
2
(ed. T. Y. T. Wu, C. J. Brokaw and C. Brennen), pp.
783
-796. New York: Plenum.
Ellington, C. P. (
1984a
). The aerodynamics of hovering insect flight. I. Quasi-steady analysis.
Phil. Trans. R. Soc. Lond. B
305
,
1
-15.
Ellington, C. P. (
1984b
). The aerodynamics of hovering insect flight. II. Morphological parameters.
Phil. Trans. R. Soc. Lond. B
305
,
17
-40.
Ellington, C. P. (
1984c
). 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.
Ennos, A. R. (
1989
). The kinematics and aerodynamics of the free flight of some Diptera.
J. Exp. Biol.
142
,
49
-85.
Fry, S. N., Sayaman, R. and Dickinson, M. H.(
2003
). The aerodynamics of free-flight maneuvers in Drosophila.
Science
300
,
495
-498.
Hamdani, H. and Sun, M. (
2000
). Aerodynamic forces and flow structures of an airfoil in some unsteady motions at small Reynolds number.
Acta Mech.
145
,
173
-187.
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.
Lehmann, F. O. and Dickinson, M. H. (
1998
). The control of wing kinematics and flight forces in fruit flies(Drosophila spp.).
J. Exp. Biol.
201
,
384
-401.
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.
Ramamurti, R. and Sandberg, W. C. (
2002
). A three-dimensional computational study of the aerodynamic mechanisms of insect flight.
J. Exp. Biol.
205
,
1507
-1518.
Sane, S. P. and Dickinson, M. H. (
2001
). The control of flight force by a flapping wing: lift and drag production.
J. Exp. Biol.
204
,
2607
-2626.
Sane, S. P. and Dickinson, M. H. (
2002
). The aerodynamic effects of wing rotation and a revised quasi-steady model of flapping flight.
J. Exp. Biol.
205
,
1087
-1096.
Spedding, G. R. (
1992
). The aerodynamics of flight. In
Advances in Comparative and Environmental Physiology II,Mechanics of Animal Locomotion
(ed. R. McN. Alexander), pp.
51
-111. London: Springer-Verlag.
Sun, M. and Tang, J. (
2002a
). Unsteady aerodynamic force generation by a model fruit fly wing in flapping motion.
J. Exp. Biol.
205
,
55
-70.
Sun, M. and Tang, J. (
2002b
). Lift and power requirements of hovering flight in Drosophila.
J. Exp. Biol.
205
,
2413
-2427.
Sun, M. and Wu, J. H. (
2003
). Aerodynamic force generation and power requirements in forward flight in a fruit fly with modeled wing motion.
J. Exp. Biol.
206
,
3065
-3083.
Usherwood, J. R. and Ellington, C. P. (
2002a
). The aerodynamics of revolving wings. I. Model hawkmoth wings.
J. Exp. Biol.
205
,
1547
-1564.
Usherwood, J. R. and Ellington, C. P. (
2002b
). The aerodynamics of revolving wings. II. Propeller force coefficients from mayfly to quail.
J. Exp. Biol.
205
,
1565
-1576.
van den Berg, C. and Ellington, C. P. (
1997a
). The vortex wake of a `hovering' model hawkmoth.
Phil. Trans. R. Soc. Lond. B
352
,
317
-328.
van den Berg, C. and Ellington, C. P. (
1997b
). The three-dimensional leading-edge vortex of a `hovering' model hawkmoth.
Phil. Trans. R. Soc. Lond. B
352
,
329
-340.
Vogel, S. (
1967
). Flight in Drosophila. III. Aerodynamic characteristics of fly wings and wing models.
J. Exp. Biol.
46
,
431
-443.
Wang, Z. J. (
2000
). Two dimensional mechanism for insect hovering.
Physical. Rev. Lett.
85
,
2216
-2219.
Weis-Fogh, T. (
1973
). Quick estimates of flight fitness in hovering animals, including novel mechanism for lift production.
J. Exp. Biol.
59
,
169
-230.