Wing motion measurement and aerodynamics of hovering true hoverflies

Hovering is an important type of flight of insects. Most hovering insects flap their wings in a horizontal plane, called `normal hovering' (Weis-Fogh, 1973). But some of the best hoverers (e.g. true hoverflies and dragonflies who can remain motionless at a point in the air for a long time) hover with an inclined stroke plane.

In recent years, much work has been done on aerodynamics, energetics and dynamic flight stability of insect flight and considerable progress has been achieved (e.g. Ellington et al., 1996; Liu et al., 1998; Dickinson et al., 1999; Wang, 2000; Sane and Dickinson, 2001; Sun and Tang, 2002a; Taylor and Thomas, 2003; Wang et al., 2004; Sun and Xiong, 2005; Dickson et al., 2008; Wu et al., 2009; Bergou et al., 2010; Walker et al., 2010). Most of the previous studies, however, have focused on, or are related to, normal hovering, and studies on hovering with an inclined stroke plane are very few. One reason for this is that there is a lack of wing motion data on this type of hovering.

There have been some studies on measuring the wing motion in insects hovering with an inclined stroke plane. Norberg measured the wing-tip kinematics of a freely hovering dragonfly using one high-speed camera, and obtained the stroke plane angle, stroke frequency, stroke angle and elevation angle (Norberg, 1975). Wakeling and Ellington made similar measurement for dragonflies and damselflies in near-hovering and forward flight (Wakeling and Ellington, 1997). Ellington measured the wing-tip kinematics of a true hoverfly in free hovering flight using one camera, and also obtained data on stroke plane angle, stroke frequency, stroke angle and elevation angle (Ellington, 1984b). In these works (Norberg, 1975; Wakeling and Ellington, 1997; Ellington, 1984b), because only one camera was used, the continuous time variation of wing orientation (wing angle of attack, wing rotation rate at stroke reversal, etc.) could not be obtained. Wing angle of attack and its rate are important parameters for determining aerodynamic force, and without data on these parameters, wing motion could not be described. Recently, Wang and Russell obtained the continuous time variation of wing angle of attack, in addition to the stroke angle and elevation angle, of hovering dragonflies using three high-speed cameras (Wang and Russell, 2007). However, because of difficulty in handling the insects, only the case of tethered flight was considered.

Because, as aforementioned, some excellent hoverers hover with an inclined stroke plane, it is of great interest to obtain detailed free-flight wing kinematic data and study the aerodynamic mechanisms, power requirements and other problems related to this type of hovering. In the present study, we measured the time course of the three-dimensional (3-D) wing motion of hovering true hoverflies using three orthogonally aligned high-speed cameras and the required morphological data. On the basis of the measured data, we used computational fluid dynamics (CFD) to compute the aerodynamic forces of the flapping wings and the power required for hovering, and to compute the aerodynamic derivatives and analyse the dynamic properties of the flight system. Comparison between the computed results and the force balance condition (vertical force being equal to the weight) can provide a test of the computational model; analyzing the time courses of the wing motion, aerodynamic forces and flow fields can provide insights into the aerodynamic force production mechanism; comparing the power requirement of the inclined stroke-plane hovering with that of normal hovering could answer the interesting question of whether power required for these two types of hovering is very different; and comparing the dynamic flight stability properties of the inclined stroke-plane hovering with that of normal hovering might provide insight into the maneuverability associated with the two types of hovering.

We netted hoverflies Episyrphus baltealus (De Geer, 1776) in a suburb of Beijing in June 2010. All experiments were conducted on the same day of capture. Only the most vigorous individuals were selected for the experiments.

We filmed the hovering flight of the hoverflies in an enclosed flight chamber using three orthogonally aligned synchronized high-speed cameras (MotionXtra HG-LE, Redlake MASD, Inc., San Diego, CA, USA; 5000 frames s^–1, shutter speed 50 μs, resolution 512×320 pixels) mounted on an optical table. The flight chamber, a cube of 15×15×15 cm³, was built from translucent glass. We backlit each camera view using densely packed arrays of light emitting diode (LED) covered with diffusion paper. We used LED arrays as the light source because they produced much less heat than cine lights at the elevated light levels required for high-speed filming. We manually triggered the synchronized cameras when the insect was observed to hover in the approximately 5×5×5 cm³ cubic zone at the central region of the chamber, which represented the intersecting field of views of the three cameras.

After the flight events were recorded, a method was required to extract the 3-D body and wing kinematics from the filmed data. The method we used here was the same as that used by Fry and others (Fry et al., 2005; Liu and Sun, 2008) [for recent development of the measuring methods, see Belhaoua et al. (Belhaoua et al., 2009), Ristroph et al. (Ristroph et al., 2009) and Walker et al. (Walker et al., 2009)]. The method is described in detail in the Appendix and is outlined here. The body and wings were represented by models (see Appendix), e.g. the model of a wing was the outline of the wing obtained by scanning the cut-off wing (Fig. 1A). We developed an interactive graphic user interface using MATLAB (v. 7.1, The Mathworks, Inc., Natick, MA, USA) to extract the 3-D body and wing positions from the frames recorded by the three cameras (an example of the frames is shown in Fig. 1B). The positions and orientations of the models of the body and wings were adjusted until the best overlap between a model image and the displayed frame was achieved in three views (the fitting process was done manually). At this point, the positions and orientations of these models would be taken as the positions and orientations of the body and the wings.

Errors in the method came from several sources. One was the deformation of the flapping wings: the wings were not flat-plate wings as they were modeled in the method (as described in the Appendix and mentioned above, the model of a wing was the outline of the wing obtained by scanning the cut-off wing, which was necessarily of a flat plate). Other error sources included errors due to camera model inaccuracy (a pinhole model was used to characterize the camera), camera calibration, stereo rig calibration, and discretization. It was shown that errors due to camera model inaccuracy, camera calibration and stereo rig calibration were small (see Appendix), and the primary errors of the method were errors due to wing deformation and discretization. We estimated the wing deformation and discretization errors as a whole by applying the method to a computer-generated virtual insect, which had similar wing motion as the hoverflies (see Appendix). Based on the observation of the flight of many insects (Ellington, 1984b; Ennos, 1989a; Walker et al., 2009), the wing deformation of the virtual insect was assumed to have a 15 deg twist and 6% camber during the translation phase of the downstroke or upstroke and the twist and camber increased to 25 deg and 10%, respectively, at stroke reversal. Analysis showed that errors in positional angle and elevation angle of the wing were within 3 deg and errors in pitch angle (or angle of attack) were within 4 deg (see Appendix).

The present method of measuring the morphological parameters follows, for the most part, that given by Ellington, whose paper can be consulted for a more detailed description of the method (Ellington, 1984a).

The insects were killed with ethyl acetate vapor after filming. The total mass (m) was measured to an accuracy of ±0.01 mg. The wings were then cut from the body and the mass of the wingless body was measured. The wing mass (m_wg) was determined from the difference between the total mass and the mass of wingless body.

Immediately after the wings were cut from the body, the shape of one of them was scanned using a scanner (HP scanjet 4370; resolution 3600×3600 d.p.i.). A sample of the scanned picture of a wing is shown in Fig. 1A. Using the scanned picture, wing length (R, the distance between the wing base and the wing tip) and local wing chord length were measured to an accuracy greater than ±0.5%. Parameters including wing area, mean chord length and radius of second moment of wing area were computed using the measured wing shape.

The wingless body was scanned from two perpendicular directions (the dorsoventral and lateral views; Fig. 2). Following Ellington (Ellington, 1984a; Ellington, 1984b), the cross-section of the body was taken as an ellipse and a uniform density was assumed for the body. With these assumptions and the measured body shape, the center of mass of body could be estimated. Body length (l_b) and distance between the wing roots (l_r) were measured from the dorsal views; distance between the wing-base axis and the center of mass (l₁) and distance between the wing-base axis and the long axis of the body (h₁) were measured from the lateral view (using these data, moments and products of inertia about the center of mass were also estimated).

The aerodynamic forces and moments were computed using the CFD method. Under hovering flight conditions, Aono and others showed that interaction between wing and body was negligibly small: the aerodynamic force in the case with the body–wing interaction was less than 2% different from that without body–wing interaction (Aono et al., 2008; Yu and Sun, 2009). Although the left and right wings might interact via a `clap and fling' mechanism, this mechanism was irrelevant in the present study because of small stroke amplitude. Therefore, in the present CFD model, the body was neglected and only the flows around one wing were computed (the aerodynamic forces produced by the other wing were derived from the results of the computed wing). Recently, Walker et al. measured the deformation of the wings of freely hovering hoverflies (Walker et al., 2010). Using Walker et al.'s data, Du and Sun investigated the effect of wing deformation on aerodynamic forces in hovering hoverflies and showed that as a first approximation, the deformable wing could be modeled by a rigid flat-plate wing with its angle of attack being equal to the local angle of attack at the radius of second moment of wing area (Du and Sun, 2010). Thus in the CFD model, we further assumed that wings were rigid flat-plate wings; the planform of the wing was obtained from the measured data, and the wing section was a flat plate with a thickness of 3% of the chord length and with rounded leading and trailing edges.

The flow equations and solution method were the same as those used in Sun and Tang (Sun and Tang, 2002a) and the method was developed by Rogers and others (Rogers and Kwak, 1990; Rogers et al., 1991). In the method, the time derivatives of the momentum equations ware differenced using a second-order, three-point backward difference formula. To solve the time-discretized momentum equations for a divergence free velocity at a new time level, a pseudo-time level was introduced into the equations and a pseudo-time derivative of pressure divided by an artificial compressibility constant was introduced into the continuity equation. The resulting system of equations was iterated in pseudo-time until the pseudo-time derivative of pressure approached zero and, thus, the divergence of the velocity at the new time level approached zero. The derivatives of the viscous fluxes in the momentum equation were approximated using second-order central differences. For the derivatives of convective fluxes, upwind differencing based on the flux-difference splitting technique was used. A third-order upwind differencing was used at the interior points and a second-order upwind differencing was used at points next to boundaries. Details of this algorithm can be found in previous studies (Rogers and Kwak, 1990; Rogers et al., 1991).

The computational grids (OH type) were generated using a Poisson solver, which was based on the work of Hilgenstock (Hilgenstock, 1988). The grids will be further described in the Results and discussion section as will analysis of the convergence of solutions.

Boundary conditions were as follows. For the far-field boundary condition, at inflow boundary, the velocity components were specified as free-stream conditions (determined by flight speed), whereas pressure was extrapolated from the interior; at the outflow boundary, pressure was set equal to the free-stream static pressure and velocity was extrapolated from the interior. On the wing surface, impermeable wall and non-slip conditions were applied and the pressure was obtained through the normal component of the momentum equation written in the moving grid system.

Solving the Navier–Stokes equations yielded the fluid velocity and pressure at discretized grid points and time steps. The aerodynamic forces acting on the wing were calculated from the pressure and the viscous stress on the wing surface.

Four hoverflies hovering at the central region of the flight chamber (the zone of the intersecting field of view of the three cameras) were filmed. They were denoted as HF1, HF2, HF3 and HF4, respectively. All of the hoverflies except HF2 were female. For each of the hoverflies, film of approximately six wing strokes were digitized. Samples of the original video sequences for HF1 are presented as supplementary material Movie 1.

We determined the stroke plane in the same way as Ellington (Ellington, 1984b). In a wingbeat cycle, approximately 35 pictures were taken and the same number of points on the curve traced by the wing tip was determined. We projected all the wing-tip points of both the left and right wings in the six wingbeats onto the plane of symmetry of the insect. A linear regression line of these projected points on the plane of symmetry was then determined. The stroke plane was defined as a plane that passed the two wing roots and was parallel to the above line. This plane is tilted at an angle β to the horizontal (called the stroke plane angle). The angles determining the wing orientation are defined as follows. A line is drawn between the wing base and wing tip (see Fig. 1A and Fig. 3A). Let (X, Y, Z) be a reference frame with origin at the wing base and an X–Y plane coinciding with the stroke plane (Fig. 3). The orientation of the wing is determined by the three Euler angles: positional angle (ϕ), stroke deviation angle (θ) and pitch angle (Ψ) (see Fig. 3A), where ϕ is the angle between the projection of the line joining the wing base and wing tip onto the stroke plane and the Y-axis, θ is the angle between the line joining the wing base and wing tip and its projection onto the stroke plane, and Ψ is the angle between the local wing chord and line l (l is perpendicular to the wing span and parallel to the stroke plane). Ψ is related to the angle of attack of the wing (α) as follows: in the downstroke, α=Ψ; in the upstroke, α=180–Ψ. The measured data of these angles as functions of time for the left and the right wings of HF1 are shown in Fig. 4 (for each of the four insects, approximately six well-repeated wing strokes in which the left and right wings were moving symmetrically were captured). For a clear description of the data, we express time during a cycle as a non-dimensional parameter, t, such that t=0 at the start of a downstroke and t=1 at the end of the subsequent upstroke. Fig. 5 plots the phase-average value (mean ± s.d.) of the positional angle (averaged over six wingbeat cycles) for each hoverfly against non-dimensional time. Figs 6 and 7 give the corresponding data for the pitch angle and the deviation angle, respectively.

As seen in Figs 4, 5, 6, 7, the motion of the right wing is approximately the same as that of the left wing, as expected for hovering flight. For a given hoverfly, the positional angle shows less variation between successive wingbeats than the geometrical angle of attack and the deviation angle.

The stroke positional angle varied with time approximately as a sinusoidal function (Fig. 5). From the data, the stroke amplitude (Φ) and mean stroke angle () can be determined using the following equations: =(ϕ_max+ϕ_min)/2 and Φ=ϕ_max–ϕ_min, where ϕ_max and ϕ_min are the maximum and minimum values of ϕ, respectively (see Ellington 1984b) (Table 1). The stroke amplitude of the hoverflies ranged from 65 to 85 deg. Ellington obtained Φ for two freely hovering true hoverflies [the same species (E. baltealus) as in the present study]: one was in inclined hovering with Φ=66 deg and the other in normal hovering with Φ=95 deg (Ellington, 1984b). Recently, Walker and others obtained Φ for hoverflies of another species, Eristalis tenax, which only performed normal hovering (Walker et al., 2010; Liu and Sun, 2008); Φ ranged from approximately 70 to 115 deg (there were five individuals in Walker et al.'s experiment and three in Liu and Sun's experiment). We thus see that the stroke amplitude of hovering flies in inclined hovering is relatively small compared with that in normal hovering.

The pitch angle (Fig. 6) had very sharp variation at the stroke reversal (t=0.4–0.6, 0.9–1 and 0–0.1), but varied relatively slowly in the mid-position of the downstroke or upstroke (t=0.1–0.4 or 0.6–0.9). The angle of attack in the mid-position of the downstroke was approximately 50 deg, and that in the mid-position of the upstroke was approximately 20 deg. This is unlike hoverflies in normal hovering (Ellington, 1984b; Walker et al., 2010; Liu and Sun, 2008) and other normal-hovering insects (e.g. Ellington, 1984b), whose downstroke and upstroke angles of attack are not very different.

The stroke deviation angle (Fig. 7) was relatively small; it was higher at the beginning and the end of a downstroke or upstroke, and lower at the middle of the downstroke or upstroke, which led to a shallow U-shaped wing-tip trajectory (the mean deviation angle was approximately 6, 0, 4 and 1 deg for HF1, HF2, HF3 and HF4, respectively, and the amplitude of the deviation angle was approximately 4, 5, 5 and 6 deg for the four insects, respectively).

The wingbeat frequency (f), stroke plane angle (β) and body angle (χ) are also given in Table 1 [the body angle is the angle between the long axis of the body and the horizontal, see Ellington (Ellington, 1984b)].

The insects filmed in the present study were only in approximate hovering, i.e. some of them moved at very low velocity. The non-dimensional velocity of the body motion, denoted by advance ratio (J), was measured as the velocity of the body motion divided by the mean wing-tip speed 2ΦfR; the values of J are also included in Table 1. J was very small for HF1, HF3 and HF4, not more than 0.07. For HF2, J was 0.13 and the insect was in slow forward flight. However, its wing kinematics was not very different from that of the other insects, except for the upstroke pitch angle: for HF1, HF3 and HF4, Ψ was approximately constant in the mid-upstroke (≈0.6–0.9) and the wing started to pitch at ≈0.9, whereas for HF2, the wing started to pitch much earlier, at ≈0.7 (Fig. 6).

Morphological parameters of the insects are given in Table 2. Parameters in the table include total mass (m), wing mass (m_wg), wing length (R), wing area (S), radius of the second moment of wing area (r₂), body length (l_b), distance between the two wing roots (l_r), distance between the wing-base axis and the center of mass (l₁), distance from the anterior end of the body to the center of mass (l₂), distance between the wing-base axis and the long axis of the body (h₁), and the moments and product of inertia of the body about its center of mass, I_x,b, I_y,b, I_z,b and I_xz,b (here x is an axis along the long axis of the body, pointing forward, y is an axis pointing to the right side of the body, and z is the other axis; the origin is at the center of mass of the body).

With the measured wing kinematics and using the CFD method described above, aerodynamic forces produced by the flapping wings were computed. We used a sinusoidal function to fit the data and obtain the time course of ϕ, and the first two terms and first four terms of the Fourier series to fit the data of θ and α, respectively. Let V and H be the computed vertical and horizontal forces of a wing, respectively; let L and D be the lift and drag of a wing, respectively (wing lift is the force component perpendicular to the stroke plane and wing drag is the force component in the stroke plane and perpendicular to the wing span). The force coefficients are defined as C_V=V/0.5ρU²S, etc., where U is the mean velocity of wing at r₂ (U=2Φfr₂) and ρ is the fluid density. The Reynolds number (Re) is 320, 330, 260 and 240 for HF1, HF2, HF3 and HF4, respectively (Re is defined as Re=cU/ν, where c is mean chord length and ν is the kinematic viscosity of the air).

Before proceeding to study the flows and aerodynamic forces, we conducted a grid resolution test. Three grids were considered: 26×25×32 (in the normal direction of the wing surface, around the wing section and in the spanwise direction of the wing, respectively; first layer grid thickness was 0.004c), 52×51×65 (first layer grid thickness was 0.002c) and 100×99×130 (first layer grid thickness was 0.001c). Note that in each refinement, the grid dimension in each direction was approximately doubled. In the normal direction, the outer boundary was set at 20 chord lengths from the wing, and in the spanwise direction, the boundary was set at six chord lengths from the wing. Portions of the dense grid (100×99×130) are shown in Fig. 8. [Note that the wing tip is slightly different from that of the real wing shown in Fig. 1, the pointed tip being cut off; without the tip cut off, the wing tip would be much narrower than the middle portion of the wing and the grid near the wing tip would have very large distortion, which would make the computation less accurate; experimental (Usherwood and Ellington, 2002) and computational (Luo and Sun, 2005) studies have shown that a slight change in wing shape had little effect on the aerodynamic force coefficients.] The non-dimensional time step was 0.02 (non-dimensionalized by c/U; the effect of time step value was studied and it was found that a numerical solution effectively independent of the time step was achieved if the time step value was ≤0.02). Calculations were performed using the above grids for the hovering flight of HF1; the results are shown in Fig. 9. The first grid refinement (from grid 26×25×32 to grid 52×51×65) produced relatively large mean magnitudes of change in C_L and C_D of 0.19 and 0.18, respectively, but the second grid refinement (from grid 52×51×65 to grid 100×99×130) produced only a small change in the results (0.05 for both C_L and C_D). The ratio between the changes in C_L (0.05/0.19) and the changes in C_D (0.05/0.18) are approximately 1/4, as expected for the second-order method. Let us use the above data to give an estimate of the accuracy of the solution obtained by the largest grid (100×99×130). Suppose that the grid is further refined (doubling the grid dimension in each direction), one could expect that the changes in C_L and C_D would be approximately 0.0125 (0.05/4=0.0125). Based on the 1/4-convergence ratio, we could estimate that the solution by grid 100×99×130 has errors in C_L and C_D of 0.017 [0.0125×(4/3)=0.017]. The mean C_L and C_D values are 1.1 and 1.9, respectively. Therefore, it is estimated that, when using the 100×99×130 grid, the numerical discretization and convergence errors in the mean C_L and C_D are approximately 1%. The 100×99×130 grid was used for the present flow computations.

First, we looked at the computed mean forces and determined whether the vertical force could support the insect weight and the horizontal force was zero (time courses of the forces and flows will be examined in the next section). The computed mean vertical force and horizontal force coefficients ( and , respectively; averaged over one wingbeat cycle) are given in Table 3. The non-dimensional weight of an insect, denoted as C_G [C_G=mg/0.5ρU²(2S)] is also given in Table 3. For all four hoverflies considered, the weight balance condition was approximately met: was different from C_G by less than 15% (the computed mean vertical forces of HF1 and HF4 were 12 and 2% less than the weight, respectively; those of HF2 and HF3 were 6 and 14% greater than the weight, respectively). The horizontal force balance condition was also approximately met, although a little more poorly than the case of weight balance condition (the horizontal forces of HF1, HF2, HF3 and HF4 were 20, 6, 19 and 20% of their respective weights). Possible reasons for the discrepancies include idealization in the CFD model and the fact that the insects might not be in exactly balanced flight.

Here we look at the time courses of the aerodynamic forces and flows and examine the properties of the aerodynamic forces. Because the wing motions of the four insects were similar, only the results for one insect, HF1, are discussed in detail here. Fig. 10 shows the time courses of C_V, C_H, C_L and C_D for HF1 in one cycle.

From Fig. 10A, it is apparent that most of the vertical force, or the weight-supporting force, is produced in the downstroke (unlike in the case of normal hovering, in which the downstroke and upstroke contribute to the weight-supporting force approximately equally); this is expected for a flapping wing with highly inclined stroke plane. From the data in Fig. 10A, it is calculated that approximately 86% of the vertical force is produced in the downstroke (for hoverflies HF2, HF3 and HF4, this value is 89, 81 and 89%, respectively). The horizontal force produced in the downstroke is approximately the same as that in the upstroke, but they have opposite signs (Fig. 10B).

The vertical and horizontal forces of a wing are the results of the lift and drag of the wing. For HF1, the downstroke C_D is a little larger than C_L (Fig. 10C,D) (for hoverflies HF2, HF3 and HF4, C_L and C_D in the downstroke are about the same). Because the stroke plane is inclined, this means that the vertical force, or the weight-supporting force, produced in the downstroke is contributed by both the lift and the drag of the wing, unlike in the case of normal hovering, in which the weight-supporting force is mainly contributed by the lift of the wing. From data in Fig. 10C,D, it is calculated that approximately 51% of the vertical force is contributed by the drag (for hoverflies HF2, HF3 and HF4, this value is 36, 33 and 38%, respectively).

The corresponding flowfield is shown in Fig. 11. The leading-edge vortex does not shed in an entire downstroke (t=0–0.5), showing that the large C_L and C_D, and hence the large vertical force (C_V) in the downstroke, are mainly due to the delayed stall mechanism.

With the flows computed by the CFD method, the mechanical power of a wing (P) can be easily calculated:

Here, M_a is the aerodynamic moment about the wing root, M_i is the inertial moment and ω is the angular velocity vector of the wing. M_a is readily calculated using the force distribution obtained from the flow computation and ω is known from the measured data. The way to determine M_i has been described in our previous works (e.g. Sun and Tang, 2002b). In calculating M_i, moments and products of inertia of the wing mass are needed. They are estimated using the wing mass measured in the present study and the density distribution of a dronefly wing measured previously (Ennos, 1989b); the results are given in Table 4. With the moments and products of inertia known and wing acceleration computed from ω, M_i can be calculated.

The instantaneous non-dimensional power (C_p) and the non-dimensional aerodynamic (C_p,a) and inertial (C_p,i) power of a wing (non-dimensionalized by 0.5ρU²Sc) for HF1 are given in Fig. 12 (those for HF2, HF3 and HF4 are similar). It is interesting to note that the time course of C_p is more similar to that of C_p,i than to that of C_p,a, because the inertial power is larger than the aerodynamic power in many parts of the stroke cycle. This means that elastic energy storage could be important for the hoverflies.

Integrating P over the part of a wingbeat cycle where it is positive gives the positive work (W⁺); integrating P over the part of a wingbeat cycle where it is negative gives the negative work (W^–). The mass-specific power (P*) is determined as the mean mechanical power over a wingbeat cycle divided by the mass of the insect. In the case of zero elastic energy storage, with the added assumption that the cost for dissipating the negative work is negligible, P*=P₁*=W⁺/t_wm, where t_w is the wingbeat cycle. In the case of 100% elastic energy storage, P*=P₂*=(W⁺–|W^–|)/t_wm.

The computed mass-specific power is given in Table 5. In the case of zero elastic energy storage, P₁* is 39–46 W kg^–1 for the four hoverflies; in the case of 100% elastic energy storage, P₂* is 27–35 W kg^–1. For HF1, the largest possible effect of elastic energy storage reduces the power by approximately 29% (for HF2, HF3 and HF4, the values are 24, 24 and 37%, respectively).

As most of the wing mass is located near the axis of spanwise rotation of the wing and the pressure center is not far from the axis of spanwise rotation, it is expected that contribution by the spanwise rotation of the wing to the aerodynamic and inertial moments, and hence to the mechanical power, is small. Because the elevation angle is much smaller than the stroke angle, it is expected that contribution by the elevation of the wing to the mechanical power is also small. C_p, C_p,a and C_p,i, computed by neglecting contributions by the spanwise rotation and the elevation of the wing, are shown in Fig. 13, and values of the corresponding mass-specific power are included in Table 5 (numbers in parentheses). The contributions by the spanwise rotation and the elevation of the wing are indeed small: the mass-specific power without the contributions from spanwise rotation and elevation is only approximately 5% different from that with the contributions.

True hoverflies also perform normal hovering. It is of interest to know whether there is any significant difference in power between these two types of hovering for true hoverflies. We examined this problem in two ways. First, we compared the mass-specific power of the inclined-stroke hoverflies of the present study with that from previously published normal hovering simulation studies. Sun and Du, using Ellington's wing kinematics data (Ellington, 1984b), computed the mass-specific power of a true hoverfly in normal hovering; P₁* and P₂* were 39 and 27 W kg^–1, respectively (Sun and Du, 2003). These values are the same as those of HF1 in the present study, and a little smaller than those of HF2, HF3 and HF4. Next, we artificially changed the insects' body angle so that the stroke plane became horizontal and adjusted the downstroke and upstroke angles of attack so that the same weight-supporting force and horizontal force were produced. The corresponding power becomes 35.2–43.4 W kg^–1 for the case of zero elastic energy storage and 25.2–34.3 W kg^–1 for the case of 100% elastic energy storage (Table 6); these values are less than 10% smaller than those found for inclined stroke-plane hovering.

From the above comparisons, we see that for the true hoverflies, the power requirement for inclined stroke-plane hovering is only slightly (approximately 10%) larger than that of normal hovering.

With the present CFD model, aerodynamic derivatives of a hovering insect can be readily computed. With the aerodynamic derivatives and the measured morphological parameters (mass, moment of inertia, etc.; Table 2), the system matrix of the insect can be determined. The eigenvalues and eigenvectors of the system matrix give the stability properties of the flight, from which insight into flight maneuverability can be obtained. Here, we conduct a longitudinal stability analysis of the hovering flight of HF1 (among the four insects considered in the paper, HF1 has the largest inclined stroke-plane angle). A description of the stability analysis method can be found in previous studies (Taylor and Thomas, 2003; Sun and Xiong, 2005).

The longitudinal system matrix of a hovering insect (A) is as follows:

where g is the gravitational acceleration; X_u, Z_u, etc. are the aerodynamic derivatives [for example, X_u is the derivative of X with respect to u (X and Z are the x- and z-components of the aerodynamic force, respectively)]; M is the aerodynamic pitching moment (x and z are coordinate axes fixed on the body with the origin at the center of mass; at equilibrium flight, i.e. hovering, x is horizontal, pointing forward, and z is vertical, pointing downward); u and w are the x- and z-components of the velocity of the center of mass; and q is the pitching rate of the body. For comparison with other insects, non-dimensional quantities are used: X_u, X_w, etc. are non-dimensionalized by 0.5ρU(2S); X_q and Z_q by 0.5ρU²(2S)/f; M_u and M_w by 0.5ρU(2S)c; M_q by 0.5ρU²(2S)c/f; m by 0.5ρU(2S)/f; I_y,b by 0.5ρU²(2S)/f²; and g by Uf [we use the superscript `+' to denote a non-dimensional quantity, e.g. X⁺_u=X_u/0.5ρU²(2S)].

The computed aerodynamic derivatives of HF1 are given in Table 7. Sun and Wang conducted longitudinal stability analysis of a hoverfly in normal hovering (Sun and Wang, 2007); here, we call this hoverfly HF_norm. For comparison, the aerodynamic derivatives of HF_norm are included in Table 7. We also computed the stability derivatives with HF1 rotated to normal hovering (in the same way as done in the power analysis); the insect in this case is denoted as HF1_rot. The aerodynamic derivatives of HF1_rot are also given in Table 7.

As seen in Eqn 2, the elements of the system matrix, and hence the stability properties, are determined not only by the aerodynamic derivatives, but also by the derivatives divided by the mass or moment of inertia of the insect. The derivatives, normalized by mass or moment of inertia, i.e. X⁺_u/m⁺, Z⁺_u/m⁺, etc., are shown in Table 8. The difference in each element of the system matrix between HF1 and HF1_rot is very small; there is some difference between HF1 and HF_norm, but it is not very large. This indicates that the eigenvalues and eigenvectors, and hence the stability properties of HF1, HF1_rot and HF_norm, would not be very different.

The eigenvalues for the three cases are shown in Table 9. As anticipated, there is little difference in the eigenvalues between HF1 and HF1_rot, and only a small difference between HF1 and HF_norm (this is also true for the eigenvectors; data not shown).

The system matrix or the eigenvalues and eigenvectors represent the inherent property of a flight system: the less stable a system is, the easier it is to change its state. The above results, especially the results of HF1 and HF1_rot (the same insect in inclined stroke-plane hovering and in normal hovering), show that stability properties are almost the same for the two types of hovering. Thus, we see that the two types of hovering have little difference in maneuverability.

The method of measuring the kinematic parameters of the body and wings using the filmed data is based on stereovision-based triangulation (see Fry et al., 2003; Fry et al., 2005; Liu and Sun, 2008; Belhaoua et al., 2009; Ristroph et al., 2009; Walker et al., 2009). The general structure of a trinocular stereo vision system is shown in Fig. A1 [in the figure, (X_w, Y_w, Z_w) is the world coordinate system; (X_C1, Y_C1, Z_C1), (X_C2, Y_C2, Z_C2) and (X_C3, Y_C3, Z_C3) are the camera coordinate systems of camera 1, 2 and 3, respectively; and (u₁, v₁), (u₂, v₂) and (u₃, v₃) are the image coordinate systems of camera 1, 2 and 3, respectively]. Based on the basic principles of stereo vision, the projections of the scene point in the world coordinate system onto any image coordinate system can be calculated on the condition that the point's coordinate in the world coordinate system and projection matrices are known. The calculation procedure is as follows.

Let P be an arbitrary point, whose coordinates are X_w, Y_w and Z_w in the world frame and X_C1, Y_C1 and Z_C1 in the frame of camera 1. Let p₁ be the projective point of P on camera 1; p₁'s coordinates in the image frame of camera 1 are u₁ and v₁. The coordinate transformation between the world coordinate system and the image coordinate systems is: where M₁ is the projection matrix of camera 1 (the projection matrix is determined by camera calibration, see below). Let m1_i,j (i, j=1, 2, 3) denote the elements of M₁. Eliminating Z_C1 from Eqn A1 gives:

The coordinates of point P's projection onto the image coordinate system of camera 1 can be determined viaEqn A2. By consecutively replacing M₁ by projection matrices of camera 2 (M₂) and camera 3 (M₃) and repeating the above calculation, we can obtain all projections on the three image coordinate systems.

We use a line segment, whose length is equal to the body length of the insect, to represent the body of the insect, and the outline of the wing, obtained from scanning the wing (see Fig. 1A), to represent the wing. The line segment of the body and the outline of a wing are referred to as the models of the body and the wing, respectively. Each of them is represented as a set of points. Using the method described above, we could easily compute the projection of the point set on all the three image coordinate systems.

We put the models of the body and the wings into the world coordinate system, then change the positions and orientations of those models until the best overlap between a models' projection and the displayed frame is achieved in all three views. At this point, the positions and orientations of those models would be taken as the positions and orientations of the body and wings of the insect. Generally, several re-adjustments of each model's position and orientation are required to obtain a satisfactory overlap. The matching is done manually. We have developed a graphical user interface subroutine for MATLAB to execute the above procedures.

We used a flat glass panel with a high accuracy black-and-white checkerboard pattern printed on it to calibrate the three orthogonally aligned cameras. The calibration gave the intrinsic and extrinsic parameters of each camera (these parameters determined the transform matrix of the camera). Each image of the calibration panel contained 324 automatically identified corner points as required feature points. Fifty-five of the 90 images of the calibration panel were used for calibrating each single camera and the stereo rigs. The calibrations were carried out using a custom-written program in MATLAB.

The stereovision-based triangulation method has several types of error (Belhaoua et al., 2009). In the present case, the errors are due to camera model inaccuracy (pinhole model was used to characterize the cameras), camera calibration, stereo rig calibration and discretization (the projection of a 3-D point in the image plane is approximated to the nearest pixel). Another error source is wing deformation. Insect wings have time-dependent deformation during flapping motion. Our model of the wing is a rigid flat plate that has the same outline as the cut-off wing, and errors will arise when one tries to match the model's projection with the image of the flapping wing.

Errors due to the camera model, camera calibration and stereo rig calibration are categorized as system error. Thirty-five images of the calibration panel that had not been used in calibrating the cameras were used to estimate these errors and the results were as follows. We used the reprojected pixel error (RPE) (Walker et al., 2009) as an indicator of camera model accuracy. For all three cameras, the RPE was less than 0.16 pixels in both directions on the image plane. Adding distortion and/or skew parameters to the pinhole model did not reduce the RPE significantly. Accuracy in describing the relative orientation and position between every two cameras was 0.02 deg and 0.03 mm, respectively. In the above processes, all computations reached a sub-pixel level accuracy based on the features detected, so the results were not affected by the discretization errors, which will be considered below.

The above results show that errors due to the camera model, camera calibration and stereo rig calibration are very small. It is expected that the other two types of errors, i.e. discretization errors and errors due to wing deformation, are the primary errors of the method. Here, we estimate these two types of errors as a whole; this estimate will be taken as the error of the method.

To estimate discretization and wing deformation errors, we applied our method to a computer-generated virtual insect consisting of a rigid body and two deforming wings. The wings can rotate around their roots (being joints with three degrees of freedom) and have time-dependent twist and camber deformations. The body is unrestrained in space, having six degrees of freedom. The wing and body motions and the wing deformation of the virtual insect are as follows. The positional and elevation angles and the pitch angle at r₂ of the wing are set the same as those of HF1. Based on observations of flight in many insects (Walker et al., 2009; Ellington, 1984b; Ennos, 1989a), the wing is assumed to have a 15 deg twist and 6% camber during the translation phase of the downstroke or upstroke and the twist and camber increase to 25 deg and 10%, respectively, at stroke reversal. It is assumed that the center of mass of the body has a horizontal harmonic oscillation of 2.5 mm amplitude and the pitch angle of the body has a harmonic oscillation of 5 deg amplitude. Each wingbeat cycle contains 35 data points.

First, we used the projection matrices to project the virtual insect onto the image planes of the cameras (setting the resolution and pixel ratio the same as those of the real cameras) and obtained three nearly orthogonal projective images of the virtual insect (an example of the 3-D virtual insect and its images is given in Fig. A2). This step contains the discretization errors. Next, as was done in the real experiment, we represented the virtual insect by the body and wing models (body represented by a line segment and wing represented by its outline; each of the models was a set of points). We put the models of the body and the wings into the world coordinate system and then changed the positions and orientations of the models until the best overlap between each model's projection and the corresponding image of the virtual insect was achieved in all three views. We thus obtained the measurements of body and wing kinematics of the virtual insect.

Comparison between the measured and the imposed body and wing kinematics gives the errors. Six wingbeats were analyzed. Fig. A3 shows the differences between the measured and imposed kinematic parameters and the corresponding histograms for the body and the right wing (results for the left wing are similar). As seen in Fig. A3A, for positional angle and elevation angle of the wing, errors are within 3 deg and the residuals are nearly centered around zero, indicating that there are only small systematic deviations; for the pitch angle of the wing, errors are a little larger, within 4 deg. For the body orientation (Fig. A3B), errors in yaw and pitch angles are approximately 1 deg and error in roll angle is within 3 deg. Errors in the body position (Fig. A3C) are typically less than 0.2 mm, or 2% of body length.

 
• A

  system matrix

 
• c

  mean chord length

 
• C_D

  wing drag coefficient

 
• C_G

  weight coefficient

 
• C_H

  horizontal force coefficient

 
• C_L

  wing lift coefficient

 
• C_p

  non-dimensional power

 
• C_V

  vertical force coefficient

 
• _V

  mean vertical force coefficient

 
• D

  drag of a wing

 
• f

  stroke frequency

 
• g

  gravitational acceleration

 
• H

  horizontal force

 
• h₁

  distance from wing-root axis to long axis of the body

 
• I_x,b, I_y,b, I_z,b

  moments of inertia of the body about the center of mass of the body

 
• I_x,w, I_y,w, I_z,w

  moments of inertia of a wing about the wing root

 
• I_xz,b

  product of inertia of the body

 
• I_xz,w

  product of inertia of a wing

 
• J

  advance ratio

 
• L

  lift of a wing

 
• l₁

  distance from wing-root axis to body center of mass

 
• l₂

  distance from anterior end of body to center of mass

 
• l_b

  body length

 
• l_r

  distance between two wing roots

 
• m

  mass of an insect

 
• M

  total aerodynamic pitching moment about center of mass

 
• M_a

  aerodynamic moment of wing around wing root

 
• M_i

  inertial moment of wing around wing root

 
• M_q

  derivative of M with respect to q

 
• M_u

  derivative of M with respect to u

 
• M_w

  derivative of M with respect to w

 
• m_wg

  mass of one wing

 
• P*

  body-mass-specific power

 
• P_a

  aerodynamic power

 
• P_i

  inertial power

 
• q

  pitching angular velocity about the center of mass

 
• R

  wing length

 
• r₂

  radius of the second moment of wing area

 
• S

  area of one wing

 
• t

  time

 
• t_c

  wingbeat period

 
• 

  non-dimensional time ( =0 and 1 at the start and end of a cycle, respectively)

 
• u

  component of velocity along the x-axis

 
• U

  reference velocity (mean flapping velocity at r₂)

 
• V

  vertical force of a wing

 
• w

  component of velocity along the z-axis

 
• X

  x-component of the total aerodynamic force

 
• x, y, z

  coordinates in the body-fixed frame of reference (with origin at center of mass)

 
• X_q

  derivative of X with respect to q

 
• X_u

  derivative of X with respect to u

 
• X_w

  derivative of X with respect to w

 
• Z

  z-component of the total aerodynamic force

 
• Z_q

  derivative of Z with respect to q

 
• Z_u

  derivative of Z with respect to u

 
• Z_w

  derivative of Z with respect to w

 
• +

  non-dimensional quantity

 
• α

  angle of attack of wing

 
• β

  stroke plane angle

 
• θ

  deviation angle of wing

 
• λ

  generic notation for an eigenvalue

 
• ρ

  density of fluid

 
• ϕ

  positional angle of wing

 
• 

  mean positional angle

 
• Φ

  stroke amplitude

 
• χ

  body angle

 
• Ψ

  pitch angle of wing

 
• ω

  wing rotation velocity vector