Wing kinematics measurement and aerodynamics of hovering droneflies

To understand the aerodynamics, energetics and control of insect flight, it is necessary to know the time history of the aerodynamic forces and moments produced by the flapping wings. It is difficult, even impossible, to directly measure the forces and moments on the wings of a freely flying insect. Existing means of circumventing this limitation are to measure experimentally or to compute numerically the forces and moments on model insect wings (e.g. Dickinson et al., 1999; Usherwood and Ellington,2002a; Usherwood and Ellington, 2002b; Sun and Tang,2002a).

In order to use the experimental and computational methods to obtain the aerodynamic forces and moments and to study insect flight, measurements of wing kinematics and some morphological parameters are required. Other researchers have measured wing kinematics of many insects in free flight,using high-speed cine or video; and also measured morphological data of these insects (Ellington, 1984a; Ellington, 1984b; Dudley and Ellington, 1990; Willmott and Ellington, 1997). But since these reported studies used only one camera, the continuous time variation of wing orientation (geometrical angle of attack, wing rotation rate, etc.) could not be obtained. Recently, the time course of three-dimensional (3D) wing motion of freely flying fruit flies was measured using three orthogonally aligned, high-speed cameras(Fry et al., 2005). Measurements of 3D wing motion of other insects are of great interest, but some limitations to Fry et al.'s work meant that morphological parameters such as weight and position of center of mass could not be measured. If these data were also measured, one could use them to test the experimental and computational models (a reasonable test of the experimental and computational models is that the measured or computed vertical force approximately balances the insect weight and, in hovering flight, the horizontal force and the pitching moment about the centre of mass of the insect are approximately zero).

In the present study, we measure the time course of 3D wing motion of hovering droneflies using three orthogonally alined high-speed cameras and also measure the required morphological data, and then employ the method of computational fluid dynamics (CFD) to compute the aerodynamic forces and moments of the wings flapping with the measured kinematics. The reasons we chose droneflies are (1) droneflies might have different wing kinematics from that of fruit flies used in Fry et al.'s study, (2) they are capable of motionless hovering, and (3) they fly even in strongly lit laboratory conditions. Comparison between computed results and the conditions of vertical force being equal to weight, horizontal force being equal to zero, and moment about the center of mass being equal to zero, can provide a test of the computational model. Analyzing the time courses of the wing motion and aerodynamic forces and moments could provide insights into how the aerodynamic forces and moments are produced and how the flight balance is achieved. Using the computed results, the power requirements can be readily estimated and the aerodynamic efficiency of flight examined.

We netted droneflies of the species Eristalis tenax Linnaeus in a suburb of Beijing in September 2006. All experiments were conducted in the same day of capture. Only the most vigorous individuals were selected for filming.

We filmed the hovering flight of the droneflies in an enclosed flight chamber using three orthogonally aligned synchronized high-speed cameras(MotionXtra HG-LE, 5000 frames s^–1, shutter speed 50 μs,resolution 512×320 pixels). 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 diodes(LED) covered with diffusion paper. The reason we used LED arrays as the light source was that they produced much less heat than cine lights at the elevated light levels required for high-speed filming. The filming was started only when the insect hovered in the cubic zone of approximately 5×5×5cm³ at the central region of the chamber, which represented the intersecting field of views of the the three cameras.

Using a flat calibration panel to calibrate the cameras, one can obtain the transform matrices between the world coordinate system and three image coordinate systems, which are called projection matrices. 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, provided that the point's coordinate in the world coordinate system and projection matrices are known. The general structure of a trinocular stereo vision system is shown in Fig. 1.

We take a line segment, whose length is equal to the body length of the insect, to represent the body of the insect, and used the outline of the wing obtained from scanned picture of the wing (see Fig. 2A) to represent the wing. The line segment of the body and the outline of a wing are designated as the models of the body and the wing, respectively. They can be represented as a set of points. Using the method described above, we can easily computed the projection of the point set on all the three image coordinate systems.

We developed a toolbox for Matlab (The MathWorks, Inc., Natick, MA, USA) to extract the position and attitude of the body and of the wings from images captured by all three cameras. We put the models of the body and the wings into the world coordinate system, then changed the positions and attitudes of those models until the best overlap between a model's projection and displayed frame was achieved in all three views. At this point, the positions and attitudes of those models would be taken as the positions and attitudes of the body and wings of the insect. Generally, several readjustments of each model's position and attitude were required to obtain a satisfactory overlap.

The present method of measuring the morphological parameters follows, for the most part, the detailed description of the method given by Ellington(Ellington, 1984a).

The insect was 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 measured. The wing mass m_wg was determined from the difference between the total mass and the mass of wingless body.

Immediately after cutting the wings 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. 2A. Using the scanned picture, wing length R (the distance between the wing base and the wing tip) (see Ellington,1984a) and local wing chord length were measured to an accuracy better than ±0.5%. Parameters including wing area, mean chord length,radius of second moment of wing area, etc., were computed using the measured wing shape.

The wingless body (with legs removed) was scanned from two perpendicular directions (the dorsoventral and lateral views; Fig. 3). Following Ellington(Ellington, 1984a), the cross section of the body was taken as 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 determined. 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(the pitch moment of inertia about the center of mass could also be determined).

The aerodynamic forces and moments were computed using the CFD method. On the basis of studies on wing–wing interactions(Lehmann et al., 2005; Sun and Yu, 2006), aerodynamic interactions between the left and right wings could be neglected. During hovering flight, the body did not move, and it was assumed that the aerodynamic interaction between the body and the wings was negligible. Therefore in the present CFD model, the body was neglected and only the flows around one wing were computed (the aerodynamic force and moment produced by the other wing were derived from the results of the computed wing). The wing planform of the insect was obtained from the present measured data. The wing section was assumed to be a flat plate with rounded leading and trailing edges, the thickness of which was 3% of the mean chord length of the wing.

The flow equations and the solution method used were the same as those described in Sun and Tang (Sun and Tang,2002a). The computational grid had dimensions 100×99×105 in the normal direction, around the wing section and in the spanwise direction, respectively (portions of the grid used for one of the insects in the present study are shown in Fig. 4). The normal grid spacing at the wall was 0.0015c (where c is the mean chord length of wing). The outer boundary was set at 20 chord lengths from the wing. The non-dimensional time step was 0.02 (non-dimensionalized by c/U, where U is the mean speed of wing at the radius of second moment of wing area; there are about 400 time steps in a stroke cycle). A detailed study of the numerical variables such as grid size,domain size, time step, etc., was conducted and it was shown that the above values for the numerical variables were appropriate for the calculations.

Three droneflies hovering at the central region of the flight chamber (the zone of the interesting field of view of the three cameras) were filmed. They were denoted as DF1, DF2 and DF3, respectively. For each of the droneflies,film of around 12 wing strokes were digitized. Samples of the original video sequences for DF1 are presented as Movie 1 in supplementary material.

Morphological parameters of insects filmed in free hovering are given in Table 1. Parameters in the table include the total mass of an insect (m), the mass of a wing(m_wg), the wing length (R), mean chord length of wing (c), area of a wing (S), radius of second moment of wing area (r₂), body length (l_b) and distance between the two wing roots (l_r), distance between the wing-base axis and the center of mass (l₁), distance between the wing-base axis and the long axis of the body(h₁), and the pitch moment of inertia about the center of mass (I_y).

We determine the stroke plane in the same way as that of Ellington(Ellington, 1984b). The stroke angle (ϕ), the stroke deviation angle (θ) and the angle of attack(α) of a wing are defined as in Fig. 5. The measured data of these angles as functions of time for the left and the right wings of DF1 are shown in Fig. 6 (for the three insects in hovering, about seven well-repeated wing strokes in which the left and right wings moving symmetrically are captured). As seen in Fig. 6, the motion of the right wing is approximately the same as that of the left wing, as expected for hovering flight. The stroke positional angle varies with time approximately as a sinusoidal function. The angle of attack does not vary significantly in the mid-position of the down- or upstroke, but varies sharply during the stroke reversal. The stroke deviation angle is relatively small; it is higher at the begining and the end of a downstroke or upstroke, and lower at the middle of the down- or upstroke.

Results for the other two insects are similar. For all the three droneflies, parameters specifying the wing motion are determined from the data. They include:

The filmed insects are only in approximate hovering flight; i.e. some of them move at very small velocity, and the non-dimensional velocity of motion,denoted by advance ratio (J) was measured (J is the velocity of motion divided by the mean wing-tip speed 2ΦnR). The values of J are also included in Table 2 and are very small for the three insects.

With measured wing kinematics and using the method described in Materials and methods, aerodynamic forces and moments produced by the flapping wings were computed. Since the wing motions of the three insects are similar, only the computed results for one insect, DF1, are discussed in detail here.

The downstroke produces a little more wing lift and drag than the upstroke(see Fig. 8D,E), because the angle of attack of the downstroke is a little larger than that of the upstroke(α_d=33.8° and α_u=33.0°). The vertical force is approximately the same as the wing lift (comparing Fig. 8A with D), because the stroke plane is almost horizontal (β is only 4°). From the data in Fig. 8A, the azimuthal position of the line of action of the vertical force of a downstroke or an upstroke can be calculated. The calculation shows that the line of action is near the position of the mean stroke angle, about 3° before and 5° after the position of the mean stroke angle for the downstroke and the upstroke,respectively.

As seen from the time courses of force coefficients(Fig. 8), the major part of the mean vertical force, or the mean lift, come from the mid-portions of the down-and upstrokes (the large force peaks at t̂=0.12–0.38 and t̂=0.62–0.88) and the contribution from wing rotation around stroke reversal is relatively small(the relatively small force peaks at t̂=0.38–0.5 and t̂=0.88–1). As discussed by Sun and Tang, the large C_L in the mid-portion of a down- or upstroke (during which period the wing is in pure translational motion) is due to the delayed stall mechanism, and the relatively small C_L peak near the end of the downstroke or upstroke is due to the fast pitching-up rotation of wing(Sun and Tang, 2002a). Using the data in Fig. 8, it is estimated that around 60% of the mean vertical force is contributed by the pure translational motion, or is due to the delayed stall mechanism.

The flow-field data can provide further evidence for the above statement. The contours of the non-dimensional spanwise component of vorticity at mid-span location are given in Fig. 9. The leading-edge vortex does not shed in an entire down- or upstroke, showing that the large C_L and C_D in the mid-portion of the half-stroke are due to the delayed stall mechanism.

The instantaneous non-dimensional power (C_p),non-dimensional aerodynamic (C_p,a) and inertial power(C_p,i) of a wing (non-dimensionalized by 0.5ρU²Sc) for DF1 are given in Fig. 10 (those for DF2 and DF3 are similar). It is interesting to note that the time course of C_p is similar to that of C_p,i; this is because the inertial power is large than the aerodynamic power. This means that elastic energy storage could be important for the droneflies.

From the measured data of wing kinematics(Table 2), we have the following main observations. The wing beat frequency (n) of the three hovering droneflies ranges from 142 Hz to 209 Hz. The stroke amplitude (Φ)is around 110°. The downstroke angle of attack (α_d) is around 35° [the upstroke angle of attack (α_u) varies with the stroke plane angle; see below]. The stroke plane angle (β) ranges from 4° to 16.4°. When β is small, α_u is about the same as α_d (e.g. for DF1, is 4°, α_uis 33.0°, approximately the same as α_d=33.8°); whenβ is relatively large (for DF2, is about 16°), α_u is smaller than α_d by about 10°. Wing rotation time for spination is a little longer than that of pronation (e.g. for DF1, the former is 32% of the wing stroke period and the latter is 26%). Wing rotation is approximately symmetrical. The mean stroke angle (⁠

Using one cine camera, Ellington measured the positional angle and deviation angle of many insects, including droneflies(Ellington, 1984b). His data show that the ϕ varies with time approximately as a simple harmonic function, and that θ is small and the wing-tip trajectory is a shallow U-shaped curve. Our data of ϕ, θ and the wing-tip trajectory agree with Ellington's [comparing Fig. 7A and Fig. 11A of the present paper with fig. 13of Ellington (Ellington,1984b)].

Because only one camera was used, Ellington could not measure the angle of attack. The present study has provided the time course of the angle of attack,giving data on angles of attack during downstroke and upstroke translations and rotation duration timing.

In a previous study, Fry et al. measured the wing kinematics of fruit flies using the same method as we use here. It is of great interest to compare the kinematics of the droneflies measured here with that of the fruit flies. For comparison, we replotted the data of the fruit flies in Fig. 12 [the data are taken from fig. 2 of Fry et al.(Fry et al., 2005); note that in Fry et al.'s paper, α is zero when the wing surface is perpendicular to the stroke plane, while in the present study α is 90° when wing is in this position, and we have converted the α values of Fry et al. to fit our definition]; data from DF1 are also given in the figure.

For both insects, ϕ varies with time approximately according to a simple harmonic function, except that for the fruit fly the duration of downstroke is slightly longer than that of the upstroke (comparing Fig. 12A with Fig. 12D). However, large differences in α and θ exist between the two insects.

For the fruit fly, after the supination (downstroke/upstroke) rotation and at the beginning of the following upstroke(Fig. 12B, t̂≈0.6), the α value is nearly 170°, i.e. the wing chord is almost parallel to the horizontal plane, and then α decreases to about 130° near the end of the upstroke; αvaries throughout the half-stroke. But for the dronefly, after the supination,α approximately keeps a constant value (around 145°) in the upstroke(Fig. 12E).

For the fruit fly, stroke deviation angle is larger, peak-to-peak amplitude is around 30°, whilst for the dronefly, this value is around 10°. As a result, the wing-tip trajectory of the fruit fly is a deep U-shaped curve and that of the dronefly a very shallow one (compare Fig. 11A, B). Moreover,θ of the fruit fly has a fast and large decrease at the beginning of the upstroke (Fig. 12C, at t̂≈0.6), shifting the minimum point of the U-shaped curve of the upstroke forward (see Fig. 11B); but for dronefly,the minimum point of the U-shaped curve is approximately at the middle (see Fig. 11A).

As discussed below, the differences in α and θ between the two insects have significant effects on how the flight is balanced.

In their study of hovering fruit fly flight, Fry et al. `replayed' the measured wing kinematics on a robotic model wing-pair and measured the aerodynamic forces of the wing-pair. In the present study, we used the measured wing kinematics of the droneflies in a CFD model to compute the aerodynamic forces and their corresponding moments. It is of great interest to compare the aerodynamic forces of the fruit flies and that of the droneflies.

When hovering with an approximately horizontal stroke plane, for the dronefly, the aerodynamic forces of the upstroke are a little smaller than that of the downstroke. For example, for DF1 (β=4.0°), the vertical force produced during the upstroke is approximately 11% less than that during the downstroke. By contrast, for fruit fly, the aerodynamic forces of the upstroke are larger than those of downstroke. For comparison, we have re-plotted the vertical force data of the fruit flies in Fig. 13 [using data taken from fig. 3A of Fry et al.(Fry et al., 2005)]. From the data, we estimate that the vertical force produced during the upstroke is approximately 58% larger than that during the downstroke.

For the dronefly, the vertical force that balances the insect weight is mainly produced in the mid positions of the down- and upstrokes (see Fig. 8A), during which the wing is in translational motion with constant angle of attack. Since the stroke plane is approximately horizontal and the variation of the deviation angle is small in the mid portions of the downstroke and upstroke, it can be said that the weight-supporting vertical force is mainly contributed by the lift of the wings. However, for the fruit fly, a large portion of the vertical force is produced in the early upstroke (Fig. 13). In this period of the upstroke, as discussed in Fry et al.'s paper, the wing surface is almost horizontal and the wing moves downward(Fry et al., 2005). Therefore,the vertical force produced in this period is contributed by the drag of the wings. Thus we see that when hovering with an approximately horizontal stroke plane, the dronefly mainly uses the lift principle for the weight-supporting force, whilst the fruit fly uses both lift and drag principles.

As shown above, for droneflies, the time course of mechanical power is similar to that of the inertial power (see Fig. 10), because the magnitude of inertial power is larger than the aerodynamic power. However, for the fruitflies, the inertial power is smaller than the aerodynamic power and the time course of mechanical power is similar to that of aerodynamic power[see Fig. 14 in which we have plotted the power data of the fruit flies using data from Fry et al.(Fry et al., 2005)]. As a result, droneflies could benefit from elastic energy storage much more than fruit flies: for the droneflies, the largest possible effect of elastic energy storage amounts a reduction of power by around 40%, whilst for fruit flies,the value is around 16% [calculated using data in table 1 of Fry et al.(Fry et al., 2005):(115–97)/115=0.16].

For the fruit flies, although n is a little larger than that of the droneflies [around 218 Hz (Fry et al.,2005); n for the droneflies is around 172 Hz], Rand m_wg/m are very small [R=2.39 mm(Fry et al., 2005); m_wg/m=0.12%, estimated using data of Vogel(Vogel, 1965); note that the m_wg is the mass of one wing]; these values are only about one fifth and one fourth of those of the droneflies, respectively], thus the inertial torque is small compared to the aerodynamic torque, resulting the small inertial power compared to the aerodynamic power. The large inertial power (compared to the aerodynamic power) for the droneflies can be explained by similar reasoning.

In our previous studies on flight force requirements and dynamic stability(e.g. Sun and Tang, 2002b; Sun and Wang, 2007), wing kinematic parameters for equilibrium flight, in which the forces and moments acting on the insect are balanced, must be given. Wing angles of attack(a_d and a_u) and mean stroke angle(⁠

Now that the angles of attack and the mean stroke angle have been measured in the present study (along with other wing kinematic parameters and the morphological parameters of wing and body), it is of interest to input all these parameters to the computation model and see whether or not an equilibrium flight is realized. This could provide a test of the computation model.

The non-dimensional weight of an insect, denoted as C_G,is defined as C_G=mg/0.5ρU²(2S). Using data in Tables 1 and 2, values of C_G for DF1, DF2 and DF3 are computed and given in Table 3. The results of the mean vertical force, horizontal force and pitching moment coefficients(C̄_V, C̄_H and C̄_M), calculated using the computation model, are also shown in Table 3. It is seen that for all three droneflies considered, the equilibrium flight conditions are approximately met: C̄_V is close to C_G (C̄_Vis different from C_G by less than 6%) and C̄_H and C̄_M are close to zero. These results show that the computation model is reasonably accurate.

The present measurements on droneflies show that the wing moves at a constant angle of attack during the translational phase of a half-stroke and has a shallow U-shaped wing-tip trajectory. Based on the wing kinematic pattern, analysis using the CFD model shows that the droneflies primarily use lift to produce weight-supporting force during the translational phase(via the delayed stall mechanism).

Fry et al.'s measurements on fruit flies(Fry et al., 2005) showed a very different wing kinematic pattern: the wing had a large downward plunge at the start of a half-stroke, resulting in a deep U-shaped wing-tip trajectory. Based on this wing kinematics, their experimental study using a robotic wing showed that the fruit flies rely heavily on drag mechanism during stroke reversal in producing vertical force.

These results demonstrate that insects having different wing kinematic patterns may employ different aerodynamic mechanisms for flight, and that in order to reveal the aerodynamic mechanisms an insect uses, detailed wing kinematics measurements should be conducted first. It is suggested that the studies on hovering of fruit flies (Fry et al., 2005) and droneflies (present study) are extended to other flight modes, such as forward flight and maneuver, and later on, other representative insects should be studied.

An advantage of the CFD method in the study of insect flight aerodynamics and dynamics is that when the computational code is validated for an insect in a certain flight mode (for example, hovering), it can be readily used for other flight modes and for other insects. When changing from hovering to forward flight or other flight modes, only the input boundary conditions,determined by the wing motion, need to be changed, and this can be easily accomplished. When the code is used for other insects, the wing shape, hence the computational grid, and the Reynolds number in the Navier–Stokes equations need to be changed (Wu and Sun,2004); both can be accomplished without much difficulty.

A second advantage is that the CFD method can provide any physical quantities that are needed for flow analysis. For example, aerodynamic force distribution on the wing is available, thus total aerodynamic force, moment about the center of mass of body (for flight balance study), and moment about wing root (for power calculation) can be readily obtained. Another example is that streamline patterns and vortices in the flow can be easily visualized.

The CFD model in the present study has been successfully validated for hovering droneflies. As discussed above, it can be readily used for study of other flight modes of droneflies and flight of other insects, provided that detailed wing kinematic patterns are measured. Combined with the method of 3D high-speed videography, the CFD method could be a tool of great promise in the future study of insect flight aerodynamics and dynamics.

LIST OF ABBREVIATIONS AND SYMBOLS

 
• c̄

  mean chord length

 
• C_D

  wing drag coefficient

 
• C_G

  weight coefficient

 
• C_H

  horizontal force coefficient

 
• C̄ _H

  mean horizontal force coefficient

 
• C_L

  wing lift coefficient

 
• C_M

  pitching moment coefficient

 
• C̄ _M

  mean pitching moment coefficient

 
• C_P

  non-dimensional power

 
• C_V

  vertical force coefficient

 
• C̄ _V

  mean vertical force coefficient

 
• C_w

  non-dimensional work per cycle

 
• \(C_{\mathrm{w}}^{+}\)

  non-dimensional positive work per cycle

 
• \(C_{\mathrm{w}}^{-}\)

  non-dimensional negative work per cycle

 
• CFD

  computational fluid dynamics

 
• D

  drag of a wing

 
• g

  the gravitational acceleration

 
• h₁

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

 
• H

  horizontal force

 
• I_wg

  moment of inertia of the wing about wing root

 
• I_y

  moment of inertia of the body about center of mass

 
• J

  advance ratio

 
• l_b

  body length

 
• l_r

  distance between two wing roots

 
• l₁

  distance from wing-root axis to body center of mass

 
• L

  lift of a wing

 
• LED

  light emitting diode

 
• m

  mass of an insect

 
• m_wg

  mass of one wing

 
• M_a

  aerodynamic moment of wing around wing root

 
• M_i

  inertial moment of wing around wing root

 
• M_y

  pitching moment about center of mass

 
• n

  stroke frequency

 
• P_a

  aerodynamic power

 
• P_i

  inertial power

 
• P^*

  body-mass-specific power

 
• r₂

  radius of the second moment of wing area

 
• r_2,m

  radius of the second moment of wing mass

 
• r_d

  radius of the first moment of wing drag

 
• R

  wing length

 
• S

  area of one wing

 
• t

  time

 
• t_c

  period of wing beat

 
• t̂

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

 
• U

  reference velocity (mean flapping velocity at r²)

 
• V

  vertical force of a wing

 
• α

  angle of attack

 
• αd

  constant value of α at the mid-position of the downstroke

 
• αu

  constant value of α at the mid-position of the upstroke

 
• β

  stroke plane angle

 
• Δt̂_p

  non-dimensional time duration for pronation

 
• Δt̂_p,t

  non-dimensional advance time for pronation

 
• Δt̂_s

  non-dimensional time duration for supination

 
• Δt̂_s,t

  non-dimensional advance time for spination

 
• η1

  phase angle of the first harmonic of Fourier series of deviation angle

 
• η2

  phase angle of the second harmonic of Fourier series of deviation angle

 
• θ

  deviation angle

 
• θ0

  constant value of Fourier series of deviation angle

 
• θ1

  amplitude of the first harmonic of Fourier series of deviation angle

 
• θ2

  amplitude of the second harmonic term of deviation angle

 
• ρ

  density of fluid

 
• ν

  kinematic viscosity of fluid

 
• ϕ

  positional angle

 
• \(\overline{{\phi}}\)

  mean positional angle

 
• Φ

  stroke amplitude

 
• χ

  body angle

 
• χ0

  free body angle

 
• Ω

  wing rotation velocity vector

This research was supported by the National Natural Science Foundation of China (10732030).