SUMMARY
Most hovering insects flap their wings in a horizontal plane (body having a large angle from the horizontal), called `normal hovering'. But some of the best hoverers, e.g. true hoverflies, hover with an inclined stroke plane (body being approximately horizontal). In the present paper, wing and body kinematics of four freely hovering true hoverflies were measured using threedimensional highspeed video. The measured wing kinematics was used in a Navier–Stokes solver to compute the aerodynamic forces of the insects. The stroke amplitude of the hoverflies was relatively small, ranging from 65 to 85 deg, compared with that of normal hovering. The angle of attack in the downstroke (∼50 deg) was much larger that in the upstroke (∼20 deg), unlike normalhovering insects, whose downstroke and upstroke angles of attack are not very different. The major part of the weightsupporting force (approximately 86%) was produced in the downstroke and it was contributed by both the lift and the drag of the wing, unlike the normalhovering case in which the weightsupporting force is approximately equally contributed by the two halfstrokes and the lift principle is mainly used to produce the force. The massspecific power was 38.59–46.3 and 27.5–35.4 W kg^{–1} in the cases of 0 and 100% elastic energy storage, respectively. Comparisons with previously published results of a normalhovering true hoverfly and with results obtained by artificially making the insects' stroke planes horizontal show that for the true hoverflies, the power requirement for inclined strokeplane hover is only a little (<10%) larger than that of normal hovering.
INTRODUCTION
Hovering is an important type of flight of insects. Most hovering insects flap their wings in a horizontal plane, called `normal hovering' (WeisFogh, 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 wingtip kinematics of a freely hovering dragonfly using one highspeed 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 nearhovering and forward flight (Wakeling and Ellington, 1997). Ellington measured the wingtip 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 highspeed 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 freeflight 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 threedimensional (3D) wing motion of hovering true hoverflies using three orthogonally aligned highspeed 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 strokeplane 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 strokeplane hovering with that of normal hovering might provide insight into the maneuverability associated with the two types of hovering.
MATERIALS AND METHODS
Animals and 3D highspeed filming
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 highspeed cameras (MotionXtra HGLE, 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^{3}, 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 highspeed filming. We manually triggered the synchronized cameras when the insect was observed to hover in the approximately 5×5×5 cm^{3} cubic zone at the central region of the chamber, which represented the intersecting field of views of the three cameras.
Measurement of wing and body kinematics
After the flight events were recorded, a method was required to extract the 3D 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 cutoff wing (Fig. 1A). We developed an interactive graphic user interface using MATLAB (v. 7.1, The Mathworks, Inc., Natick, MA, USA) to extract the 3D 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 flatplate 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 cutoff 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 computergenerated 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).
Measurement of morphological parameters
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 crosssection 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 wingbase axis and the center of mass (l_{1}) and distance between the wingbase axis and the long axis of the body (h_{1}) were measured from the lateral view (using these data, moments and products of inertia about the center of mass were also estimated).
Computation of aerodynamic forces and power requirements
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 flatplate 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 flatplate 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 secondorder, threepoint backward difference formula. To solve the timediscretized momentum equations for a divergence free velocity at a new time level, a pseudotime level was introduced into the equations and a pseudotime derivative of pressure divided by an artificial compressibility constant was introduced into the continuity equation. The resulting system of equations was iterated in pseudotime until the pseudotime 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 secondorder central differences. For the derivatives of convective fluxes, upwind differencing based on the fluxdifference splitting technique was used. A thirdorder upwind differencing was used at the interior points and a secondorder 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 farfield boundary condition, at inflow boundary, the velocity components were specified as freestream conditions (determined by flight speed), whereas pressure was extrapolated from the interior; at the outflow boundary, pressure was set equal to the freestream static pressure and velocity was extrapolated from the interior. On the wing surface, impermeable wall and nonslip 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.
RESULTS AND DISCUSSION
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.
Wing kinematics
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 wingtip 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 Yaxis, θ 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 wellrepeated 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 nondimensional 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 phaseaverage value (mean ± s.d.) of the positional angle (averaged over six wingbeat cycles) for each hoverfly against nondimensional 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 midposition of the downstroke or upstroke (t=0.1–0.4 or 0.6–0.9). The angle of attack in the midposition of the downstroke was approximately 50 deg, and that in the midposition 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 normalhovering 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 Ushaped wingtip 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 nondimensional velocity of the body motion, denoted by advance ratio (J), was measured as the velocity of the body motion divided by the mean wingtip 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 midupstroke (≈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
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_{2}), body length (l_{b}), distance between the two wing roots (l_{r}), distance between the wingbase axis and the center of mass (l_{1}), distance from the anterior end of the body to the center of mass (l_{2}), distance between the wingbase axis and the long axis of the body (h_{1}), 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).
Computed aerodynamic forces
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^{2}S, etc., where U is the mean velocity of wing at r_{2} (U=2Φfr_{2}) 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 nondimensional time step was 0.02 (nondimensionalized 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 secondorder 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/4convergence 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 nondimensional weight of an insect, denoted as C_{G} [C_{G}=mg/0.5ρU^{2}(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.
Time courses of the aerodynamic forces and flows
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 weightsupporting force, is produced in the downstroke (unlike in the case of normal hovering, in which the downstroke and upstroke contribute to the weightsupporting 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 weightsupporting 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 weightsupporting 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 leadingedge 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.
Flight power
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.
(1)The instantaneous nondimensional power (C_{p}) and the nondimensional aerodynamic (C_{p,a}) and inertial (C_{p,i}) power of a wing (nondimensionalized by 0.5ρU^{2}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 massspecific 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_{1}*=W^{+}/t_{w}m, where t_{w} is the wingbeat cycle. In the case of 100% elastic energy storage, P*=P_{2}*=(W^{+}–W^{–})/t_{w}m.
The computed massspecific power is given in Table 5. In the case of zero elastic energy storage, P_{1}* is 39–46 W kg^{–1} for the four hoverflies; in the case of 100% elastic energy storage, P_{2}* 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 massspecific 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 massspecific 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 massspecific power of the inclinedstroke 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 massspecific power of a true hoverfly in normal hovering; P_{1}* and P_{2}* 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 weightsupporting 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 strokeplane hovering.
From the above comparisons, we see that for the true hoverflies, the power requirement for inclined strokeplane hovering is only slightly (approximately 10%) larger than that of normal hovering.
Aerodynamic derivatives and flight stability analysis
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 strokeplane 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 zcomponents 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 zcomponents of the velocity of the center of mass; and q is the pitching rate of the body. For comparison with other insects, nondimensional quantities are used: X_{u}, X_{w}, etc. are nondimensionalized by 0.5ρU(2S); X_{q} and Z_{q} by 0.5ρU^{2}(2S)/f; M_{u} and M_{w} by 0.5ρU(2S)c; M_{q} by 0.5ρU^{2}(2S)c/f; m by 0.5ρU(2S)/f; I_{y}_{,b} by 0.5ρU^{2}(2S)/f^{2}; and g by Uf [we use the superscript `+' to denote a nondimensional quantity, e.g. X ^{+}_{u}=X_{u}/0.5ρU^{2}(2S)].
(2)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 strokeplane 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.
APPENDIX
Measurement method and error assessment
Measuring wing and body kinematics
The method of measuring the kinematic parameters of the body and wings using the filmed data is based on stereovisionbased 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_{1}, v_{1}), (u_{2}, v_{2}) and (u_{3}, v_{3}) 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.
(A1)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_{1} be the projective point of P on camera 1; p_{1}'s coordinates in the image frame of camera 1 are u_{1} and v_{1}. The coordinate transformation between the world coordinate system and the image coordinate systems is: where M_{1} 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_{1}. Eliminating Z_{C1} from Eqn A1 gives:
(A2)The coordinates of point P's projection onto the image coordinate system of camera 1 can be determined via Eqn A2. By consecutively replacing M_{1} by projection matrices of camera 2 (M_{2}) and camera 3 (M_{3}) 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 readjustments 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.
Camera calibration
We used a flat glass panel with a high accuracy blackandwhite 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. Fiftyfive 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 customwritten program in MATLAB.
Error analysis
The stereovisionbased 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 3D point in the image plane is approximated to the nearest pixel). Another error source is wing deformation. Insect wings have timedependent deformation during flapping motion. Our model of the wing is a rigid flat plate that has the same outline as the cutoff 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. Thirtyfive 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 subpixel 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 computergenerated 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 timedependent 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_{2} 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 3D 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.
FOOTNOTES

Supplementary material available online at http://jeb.biologists.org/lookup/suppl/doi:10.1242/jeb.054874//DC1

This research was supported by grants from the National Natural Science Foundation of China (10732030 and 10902011) and the 111 Project (B07009).
LIST OF SYMBOLS
 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}
 nondimensional 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_{1}
 distance from wingroot 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_{1}
 distance from wingroot axis to body center of mass
 l_{2}
 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*
 bodymassspecific power
 P_{a}
 aerodynamic power
 P_{i}
 inertial power
 q
 pitching angular velocity about the center of mass
 R
 wing length
 r_{2}
 radius of the second moment of wing area
 S
 area of one wing
 t
 time
 t_{c}
 wingbeat period
 nondimensional time ( =0 and 1 at the start and end of a cycle, respectively)
 u
 component of velocity along the xaxis
 U
 reference velocity (mean flapping velocity at r_{2})
 V
 vertical force of a wing
 w
 component of velocity along the zaxis
 X
 xcomponent of the total aerodynamic force
 x, y, z
 coordinates in the bodyfixed 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
 zcomponent 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
 +
 nondimensional 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
 © 2011.