## SUMMARY

Aerodynamic force generation and mechanical power requirements of a dragonfly (*Aeschna juncea*) in hovering flight are studied. The method of numerically solving the Navier–Stokes equations in moving overset grids is used.

When the midstroke angles of attack in the downstroke and the upstroke are set to 52° and 8°, respectively (these values are close to those observed), the mean vertical force equals the insect weight, and the mean thrust is approximately zero. There are two large vertical force peaks in one flapping cycle. One is in the first half of the cycle, which is mainly due to the hindwings in their downstroke; the other is in the second half of the cycle, which is mainly due to the forewings in their downstroke. Hovering with a large stroke plane angle (52°), the dragonfly uses drag as a major source for its weight-supporting force (approximately 65% of the total vertical force is contributed by the drag and 35% by the lift of the wings).

The vertical force coefficient of a wing is twice as large as the quasi-steady value. The interaction between the fore- and hindwings is not very strong and is detrimental to the vertical force generation. Compared with the case of a single wing in the same motion, the interaction effect reduces the vertical forces on the fore- and hindwings by 14% and 16%, respectively,of that of the corresponding single wing. The large vertical force is due to the unsteady flow effects. The mechanism of the unsteady force is that in each downstroke of the hindwing or the forewing, a new vortex ring containing downward momentum is generated, giving an upward force.

The body-mass-specific power is 37 W kg^{-1}, which is mainly contributed by the aerodynamic power.

## Introduction

Dragonflies are capable of long-time hovering, fast forward flight and quick manoeuvres. Scientists have always been fascinated by their flight. Kinematic data such as stroke amplitude, inclination of the stroke-planes,wing beat frequency and phase-relation between the fore- and hindwings were measured by taking high-speed pictures of dragonflies in free-flight (e.g. Norberg, 1975; Wakeling and Ellington, 1997b)and tethered dragonflies (e.g. Alexander,1984). Using these data in quasi-steady analyses (not including the interaction effects between forewing and hindwing), it was shown that the lift coefficient required for flight was much greater than the steady-state values measured from dragonfly wings(Wakeling and Ellington,1997a). This suggested that unsteady wing motion and/or flow interaction between the fore- and hindwings must play important roles in the flight of dragonflies (Norberg,1975; Wakeling and Ellington,1997c).

Force measurement on a tethered dragonfly was conducted by Somps and Luttges (1985). It was shown that over some part of a stroke cycle, lift force was many times larger than that measured from dragonfly wings under steady-state conditions. This clearly showed that the effect of unsteady flow and/or wing interaction was important. Flow visualization studies on flapping model dragonfly wings were conducted by Saharon and Luttges (1988, 1989), and it was shown that constructive or destructive wing/flow interactions might occur, depending on the kinematic parameters of the flapping motion. In these studies, only the total force of the fore- and hindwings was measured and, moreover, force measurements and flow visualizations were conducted in separated works.

In order to further understand the dragonfly aerodynamics, it was desirable to determine the aerodynamic force and flow structure simultaneously and also to know the force on the individual forewing and hindwing during their flapping motions. Freymuth(1990) conducted force measurement and flow visualization on an airfoil in hover modes. One of the hover modes was for hovering dragonflies. Only mean vertical force was measured. It was shown that large mean vertical force coefficient could be obtained and the force was related to a wake of vortex pairs that produced a downward jet of stream. Wang(2000) used a computational fluid dynamics (CFD) method to study the aerodynamic force and vortex wake structure of an airfoil in dragonfly hovering mode. Time variation of the aerodynamic force in each flapping cycle and the vortex shedding process was obtained. It was shown that large vertical force was produced during each downstroke and the mean vertical force was enough to support the weight of a typical dragonfly. During each downstroke, a vortex pair was created. The large vertical force was explained by the downward two-dimensional jet induced by the vortex pair.

In the works of Freymuth(1990) and Wang(2000), only a single airfoil was considered. Lan and Sun(2001c) studied two airfoils in dragonfly hovering mode using the CFD method. For comparison, they also computed the flow of a single airfoil. For the case of the single airfoil,their results of aerodynamic force and flow structure were similar to that of Freymuth's (1990) experiment and Wang's (2000) computation. For the fore and aft airfoils flapping with 180° phase difference (counter stroking), the time variation of the aerodynamic force on each airfoil was broadly similar to that of the single airfoil; the major effect of interaction between the fore and aft airfoils was that the vertical forces on both the airfoils were decreased by approximately 20% in comparison with that of the single airfoil.

The above works (Freymuth,1990; Wang, 2000; Lan and Sun, 2001c), which obtained aerodynamic force and flow structure simultaneously, were done for airfoils. It is well known that the lift on an airplane wing of large aspect ratio can be explained by a two-dimensional wing theory. But for a dragonfly wing, although its aspect ratio is relatively large, its motion is much more complex than that of an airplane wing. Three-dimensional effect should be investigated. Moreover, the effect of aerodynamic interaction between the fore- and hindwings in three-dimensional cases is unknown. The work of Lan and Sun (2001c) on two airfoils flapping with 180° phase difference showed that interaction between the two airfoils was detrimental to their aerodynamic performance. This result is opposite to the common expectation that wing interaction of a dragonfly would enhance its aerodynamic performance. It is of interest to investigate the interaction effect in the three-dimensional case.

In the present study, we extend our previous two-dimensional study(Lan and Sun, 2001c) to a three-dimensional case. As a first step, we study the case of hovering flight. For the dragonfly *Aeschna juncea* in free hovering flight, detailed kinematic data were obtained by Norberg(1975). Morphological data of the dragonfly (wing shape, wing size, wing mass distribution, weight of the insect, etc.) are also available (Norberg,1972). On the basis of these data, the flows and aerodynamic forces and the power required for producing the forces are computed and analyzed. Because of the unique feature of the motion, i.e. the forewing and the hindwing move relative to each other, the approach of solving the flow equations over moving overset grids is employed.

## Materials and methods

### The model wings and their kinematics

The fore- and hindwings of the dragonfly are approximated by two flat plates. The thickness of the model wings is 1% of *c* (where *c*is the mean chord length of the forewing). The planforms of the model wings(see Fig. 1A) are similar to those of the real wings (Norberg,1972). The two wings have the same length, but the chord length of the hindwing is larger than that of the forewing. The radius of the second moment of the forewing area is denoted by *r*_{2}(*r*_{2}=∫_{Sf}*r*^{2}d*S*_{f}/*S*_{f}, where *r* is radial distance and *S*_{f} is the area of the forewing); *r*_{2}=0.61*R* (*R* is the wing length). The flapping motions of the wings in hovering flight are sketched in Fig. 1B. The hindwing leads the forewing in phase by 180° (Norberg,1975). The azimuthal rotation of the wing about the *z*axis (see Fig. 1C) is called`translation', and the pitching (or flip) rotation of the wing near the end of a half-stroke and at the beginning of the following half-stroke is called rotation. The speed at *r*_{2} is called the translational speed.

*u*

_{t}and is given by:

*u*

^{+}

_{t}=

*u*

_{t}/

*U*(

*U*is the reference velocity); the non-dimensional timeτ=

*tU*/

*c*(

*t*is the time;

*c*is the mean chord length of the forewing, used as reference length in the present study);τ

_{c}is the non-dimensional period of the flapping cycle; andγ is the phase angle of the translation of the wing. The reference velocity is

*U*=2Φ

*nr*

_{2}, where Φ and

*n*are the stroke amplitude and stroke frequency of the forewing,respectively. Denoting the azimuth-rotational velocity as

_{d}for the downstroke and by α

_{u}for the upstroke. Around the stroke reversal, α changes with time and the angular velocity (

_{r}≤τ≤τ

_{r}+Δτ

_{r}, and the non-dimensional form

_{r}is the time at which the rotation starts; andΔτ

_{r}is the time interval over which the rotation lasts. In the time interval of Δτ

_{r}, the wing rotates fromα=α

_{d}to α=180°–α

_{u}. Therefore, when α

_{d}, α

_{u}andΔτ

_{r}are specified,

_{u}to α=α

_{d}, and the sign of the right-hand side of equation 2 should be reversed). The axis of the flip rotation is located at a distance of 1/4 chord length from the leading edge of the wing.

### The Navier–Stokes equations and solution method

*oxyz*(Fig. 1C), are as follows:

*u, v*and

*w*are three components of the non-dimensional velocity and

*p*is the non-dimensional pressure. In the non-dimensionalization,

*U, c*and

*c/U*are taken as the reference velocity, length and time, respectively.

*Re*denotes the Reynolds number and is defined as

*Re*=

*cU*/υ, where υis kinematic viscosity of the air. Equations 3, 4, 5, 6 are solved using an algorithm based on the method of artificial compressibility. The algorithm was first developed by Rogers and Kwak(1990) and Rogers et al.(1991) for a single-zone grid,and it was extended by Rogers and Pulliam(1994) to overset grids. The algorithm is outlined below.

*x,y,z*,τ) to the curvilinear coordinate system(ξ,η,ζ,τ) using a general time-dependent coordinate transformation. For a flapping wing, in order to make the transformation simple, a body-fixed coordinate system(

*o*′

*x*′

*y*′

*z*′) is also employed (Fig. 1C). In terms of the Euler angles α and φ (defined in Fig. 1C), the inertial coordinates (

*o,x,y,z*) are related to the body-fixed coordinates(

*o*′,

*x*′,

*y*′,

*z*′)through the following relationship:

_{x},ξ

_{y},ξ

_{z},ξ

_{τ}),(η

_{x},η

_{y},η

_{z},η

_{τ})and(ζ

_{x},ζ

_{y},ζ

_{z},ζ

_{τ}),which are needed in the transformed Navier–Stokes equations, can be calculated from those in the body-fixed, non-inertial coordinate system,(ξ

_{x′},ξ

_{y′},ξ

_{z′}),(η

_{x′},η

_{y′},η

_{z′})and(ζ

_{x′},ζ

_{y′},ζ

_{z′}),which need to be calculated only once. As a wing moves, the coordinate transformation functions vary with (

*x,y,z*,τ) such that the grid system moves and always fits the wing. The body-fixed non-inertial frame of reference(

*o*′,

*x*′,

*y*′,

*z*′) is used in the initial grid generation.

The time derivatives of the momentum equations are 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 is introduced into the equations, and a pseudo-time derivative of pressure divided by an artificial compressibility constant is introduced into the continuity equation. The resulting system of equations is iterated in pseudo-time until the pseudo-time derivative of pressure approaches zero; thus, the divergence of the velocity at the new time level approaches zero. The derivatives of the viscous fluxes in the momentum equation are approximated using second-order central differences. For the derivatives of convective fluxes, upwind differencing based on the flux-difference splitting technique is used. A third-order upwind differencing is used at the interior points and a second-order upwind differencing is used at points next to boundaries. Details of this algorithm can be found in Rogers and Kwak (1990) and Rogers et al. (1991). For the computation in the present work, the artificial compressibility constant is set to 100 (it has been shown that when the artificial compressibility constant varied between 10 and 300, the number of sub-iterations changes a little but the final result does not change).

With overset grids, as shown in Fig. 2, for each wing there is a body-fitted curvilinear grid, which extends a relatively short distance from the body surface; in addition, there is a background Cartesian grid, which extends to the far-field boundary of the domain (i.e. there are three grids). The solution method for a single grid is applied to each of the three grids. The wing grids capture features such as boundary layers, separated vortices and vortex/wing interactions. The background grid carries the solution to the far field. The two wing grids are overset onto the background Cartesian grid and parts of the two wing grids overlap when the two wings move close to each other. As a result of the oversetting of the grids, there are regions of holes in the wing grids and in the background grid. As the wing grids move, the holes and hole boundaries change with time. To determine the hole-fringe points, the method known as domain connectivity functions by Meakin(1993) is employed. Intergrid boundary points are the outer-boundary points of the wing grids and the hole-fringe points. Data are interpolated from one grid to another at the hole-fringe points and, similarly, at the outer-boundary points of the wing grids. In the present study, the background grid does not move and the two wing-grids move in the background grid. The wing grids are generated by using a Poisson solver that is based on the work of Hilgenstock(1988). They are of O-H type grids. The background Cartesian grid is generated algebraically. Some portions of the grids are shown in Fig. 2.

For far-field boundary conditions, at the inflow boundary, the velocity components are specified as freestream conditions while pressure is extrapolated from the interior; at the outflow boundary, pressure is set equal to the free-stream static pressure, and the velocity is extrapolated from the interior. On the wing surfaces, impermeable wall and no-slip boundary conditions are applied, and the pressure on the boundary is obtained through the normal component of the momentum equation written in the moving coordinate system. On the plane of symmetry of the dragonfly (the *XZ* plane; see Fig. 1A), flow-symmetry conditions are applied (i.e. *w* and the derivatives of *u, v*and *p* with respect to *y* are set to zero).

### Evaluation of the aerodynamic forces

*l*

_{f}and

*d*

_{f}denote the lift and drag of the forewing, respectively;

*l*

_{h}and

*d*

_{h}denote the lift and drag of the hindwing, respectively. Resolving the lift and drag into the

*Z*and

*X*directions gives the vertical force and thrust of a wing.

*L*

_{f}and

*T*

_{f}denote the vertical force and thrust of the forewing, respectively;

*L*

_{h}and

*T*

_{h}denote the vertical force and thrust of the hindwing,respectively. For the forewing:

*l*

_{f},

*d*

_{f},

*l*

_{h},

*d*

_{h},

*L*

_{f},

*T*

_{f},

*L*

_{h}and

*T*

_{h}are denoted as

*C*

_{l,f},

*C*

_{d,f},

*C*

_{l,h},

*C*

_{d,h},

*C*

_{L,f},

*C*

_{T,f},

*C*

_{L,h}and

*C*

_{T,h}, respectively. They are defined as:

*S*

_{f}and

*S*

_{h}are the areas of the fore- and hindwings, respectively. The total vertical force coefficient (

*C*

_{L}) and total thrust coefficient (

*C*

_{T}) of the fore- and hindwings are as follows:

### Data of hovering flight in Aeschna juncea

High-speed pictures of the dragonfly *Aeschna juncea* in hovering flight were taken by Norberg(1975), and the following kinematic data were obtained. For both the fore- and hingwings, the chord is almost horizontal during the downstroke (i.e. α_{d}≈β)and is close to the vertical during the upstroke; the stroke plane angle(β) is approximately 60°; the stroke frequency (*n*) is 36 Hz,the stroke amplitude (Φ) is 69°; the hindwing leads the forewing in phase by 180°. The mass of the insect (*m*) is 754 mg; forewing length is 4.74 cm;hindwing length is 4.60 cm; the mean chord lengths of the fore- and hindwings are 0.81 cm and 1.12 cm, respectively; the moment of inertial of wing-mass with respect to the fulcrum (*I*) is 4.54 g cm^{-2} for the forewing and 3.77 g cm^{-2} for the hindwing(Norberg, 1972).

On the basis of the above data, the parameters of the model wings and the wing kinematics are determined as follows. The lengths of both wings(*R*) are assumed to be 4.7 cm; the reference length (the mean chord length of the forewing) *c*=0.81 cm; the reference velocity *U*=2Φ*nr*_{2}=2.5 m s^{-1}; the Reynolds number *Re*=*Uc*/υ≈1350; the stroke periodτ _{c}=*U*/*nc*=8.58. γ is set as 180° and 0° for the fore- and hindwings, respectively. Norberg(1975) did not provide the rate of wing rotation during stroke reversal. Reavis and Luttges(1988) made measurements on similar dragonflies and it was found that maximum

^{-1}. Here,

^{-1}, giving

_{r}=3.36.

## Results and analysis

### Test of the solver

A single-grid solver based on the computational method described above was developed by Lan and Sun(2001a). It was tested by the analytical solutions of the boundary layer flow on a flat plate(Lan and Sun, 2001a) and by the measured unsteady forces on a flapping model fruit fly wing(Sun and Wu, 2003). A moving overset-grid solver for two-dimensional flow based on the above method was developed by the same authors and it was tested by comparison with the analytical solution of the starting flow around an elliptical airfoil (Lan and Sun, 2001b,c). The two-dimensional moving overset-grids solver is extended to three dimensions in the present study. The three-dimensional moving overset-grids solver is tested here in three ways. First, the flow past a starting sphere is considered, for which the approximate solution of the Navier–Stokes equations is known. Second, the code is tested by comparing with the results of the single grid. Finally, the code is tested against experimental data of a flapping model fruit fly wing by Sane and Dickinson(2001).

As a first test, it is noted that in the initial stage of the starting motion of a sphere, because the boundary layer is still very thin, the flow around the sphere can be adequately treated by potential flow theory, and the flow velocity around the sphere and the drag (added-mass force) on the sphere can be obtained analytically. The acceleration of the sphere during the initial start is a cosine function of time; after the initial start, the sphere moves at constant speed (*U*_{s}). In the numerical calculation, the Reynolds number [based on *U*_{s} and the radius (*a*) of the sphere] is set as 10^{3}. Fig. 3A shows the numerical and analytical drag coefficients (*C*_{d}) *vs*non-dimensional time (τ_{s})(*C*_{d}=*drag*/0.5ρ*U*_{s}^{2}π*a*^{2};τ _{s}=*tU*_{s}/2*a*). Betweenτ _{s}=0 and τ_{s}≈0.2, the numerical result is in very good agreement with the analytical solution. Fig. 3B shows the azimuthal velocity (*u*_{θ}) at τ_{s}=0.1 as a function of *r*/2*a* (*r* is radial distance) with fixed azimuthal angle π/2. The numerical result agrees well with the analytical solution outside the boundary layer.

In the second test, the flow around the starting sphere is computed by the single-grid code, and the results computed using the single grid and moving overset grid are compared (also in Fig. 3). They are in good agreement. For the case of the single grid,the grid is of O-O type, where the numerical coordinates (ξ,η,ζ)lie along the standard spherical coordinates. It has dimensions 100×65×129. The outer boundary is set at 30*a* from the sphere. The non-dimensional time step is 0.01. Grid sizes of 68×41×81 and 46×27×51 were also used. By comparing the results from these three grids, it was shown that the grid size of 100×65×129 was appropriate for the computation. For the case of moving overset grids, the grid system consists of two grids: one is the curvilinear grid of the sphere and the other is the background Cartesian grid. The outer boundary of the sphere grid is at 1.4*a* from the sphere surface and the out-boundary of the background grid is 30*a* from the sphere. The grid density is made similar to that of the single grid.

In the third test, the set-up of Sane and Dickinson(2001) is followed and the aerodynamic forces are computed for the flapping model fruit fly wing. The computed lift and drag are compared with the measured values in Fig. 4. In the computation, the wing grid has dimensions of 109×50×52 around the wing section, in normal direction and in spanwise direction, respectively; the outer boundary of the wing grid is approximately 2.0*c* from the wing. The background Cartesian grid has dimensions of 90×85×80 and the outer boundary is 20*c* from the wing. The non-dimensional time step is 0.02. The grid density test was conducted and it was shown that the above overset grids were appropriate for the computation. In Fig. 4A,B, the flapping amplitude is 60° and the midstroke angle of attack is 50°; in Fig. 4C,D, these quantities are 180° and 50°, respectively. The magnitude and trends with variation over time of the computed lift and drag forces are in reasonably good agreement with the measured results.

### The total vertical force and thrust; comparison with insect weight

In the calculation, the wings start the flapping motion in still air and the calculation is ended when periodicity in aerodynamic forces and flow structure is approximately reached (periodicity is reached approximately 2–3 periods after the calculation is started).

Fig. 5 shows the total vertical force and thrust coefficients in one cycle, computed by two grid systems, grid system 1 and grid system 2. In both grid systems, the outer boundary of the wing-grid was set at about 2*c* from the wing surface and that of the background grid at about 40*c* from the wings. For grid system 1, the wing grid had dimensions 29×77×45 in the normal direction, around the wing and in the spanwise direction, respectively, and the background grid had dimensions 90×72×46 in the *X*(horizontal), *Z* (vertical) and *Y* directions, respectively(Fig. 2 shows some portions of grid system 1). For grid system 2, the corresponding grid dimensions were 41×105×63 and 123×89×64. For both grid systems, grid points of the background grid concentrated in the near field of the wings where its grid density was approximately the same as that of the outer part of the wing grid. As seen in Fig. 5, there is almost no difference between the force coefficients calculated by the two grid systems. Calculations were also conducted using a larger computational domain. The domain was enlarged by adding more grid points to the outside of the background grid of grid system 2. The calculated results showed that there was no need to put the outer boundary further than that of grid system 2. It was concluded that grid system 1 was appropriate for the present study. The effect of time step value was considered and it was found that a numerical solution effectively independent of the time step was achieved if Δτ≤0.02. Therefore, Δτ=0.02 was used in the present calculations.

From Fig. 5, it is seen that there are two large *C*_{L} peaks in one cycle, one in the first half of the cycle (while the hindwing is in its downstroke) and the other in the second half of the cycle (while the forewing is in its downstroke). It should be noted that by having two large *C*_{L} peaks alternatively in the first and second halves of a cycle, the flight would be smoother. Averaging *C*_{L} (and *C*_{T}) over one cycle gives the mean vertical force coefficient (

_{d}=52° andα

_{u}=8°, respectively. These values of β,α

_{d}and α

_{u}give an approximately balanced flight and they are close to the observed values [β≈60°; during the downstroke the chord is almost horizontal (i.e.α

_{d}≈β), and during the upstroke the chord is close to vertical].

### The forces of the forewing and the hindwing

The total vertical force (or thrust) coefficient is the sum of the vertical force (or thrust) coefficient of the fore- and hindwings. Fig. 6 gives the vertical force and thrust coefficients of the fore- and hindwings. The hindwing produces a large *C*_{L,h} peak during its downstroke (the first half of the cycle) and very small *C*_{L,h} in its upstroke (the second half of the cycle); this is true for the forewing, but its downstroke is in the second half of the cycle. Comparing Fig. 6 with Fig. 5 shows that the hindwing in its downstroke is responsible for the large *C*_{L} peak in the first half of the cycle, and the forewing in its downstroke is responsible for the large *C*_{L} peak in the second half of the cycle. The contributions to the mean total vertical force by the forewing and hindwing are 42% and 58%, respectively. The vertical force on the hindwing is 1.38 times that on the forewing. Note that the area of the hindwing is 1.32 times that of the forewing. That is, the relatively large vertical force on the hindwing is mainly due to its relatively large size.

The vertical force and thrust coefficients of a wing are the results of the lift and drag coefficients of the wing. The corresponding lift and drag coefficients *C*_{l,f}, *C*_{d,f}, *C*_{l,h} and *C*_{d,h} are shown in Fig. 7. For the hindwing, *C*_{d,h} is larger than *C*_{l,h} during the downstroke of the wing, and β is large (52°). As a result, a large part of *C*_{L,h} is from *C*_{d,h}(approximately 65% of *C*_{L,h} is from *C*_{d,h} and 35% is from *C*_{l,h}). This is also true for the forewing. That is, the dragonfly uses drag as a major source for its weight-supporting force when hovering with a large stroke plane angle.

### The mechanism of the large vertical force

As shown in Fig. 6, the peak value of *C*_{L,h} is approximately 3.0 (that of *C*_{L,f} is approximately 2.6). Note that in the definition of the force coefficient, the total area of the fore- and hindwings(*S*_{f}+*S*_{h}) and the mean flapping velocity *U* are used as reference area and reference velocity, respectively. For the hindwing, if its own area and the instantaneous velocity are used as reference area and reference velocity, respectively, the peak value of the vertical force coefficient would be 3.0×[(*S*_{f}+*S*_{h})/*S*_{h}]×*U*^{2}/(π*U*/2)^{2}=2.1. Similarly, for the forewing, the peak value would be 2.6×[(*S*_{f}+*S*_{h})/*S*_{f}]×*U*^{2}/(π*U*/2)^{2}=2.4. Since the thrust coefficients *C*_{T,f} and *C*_{T,h} are small, *C*_{L,f} and *C*_{L,h} can be taken as the coefficients of the resultant aerodynamic force on the fore- and hindwings, respectively. The above shows that the peak value of resultant aerodynamic force coefficient for the forewing or hindwing is 2.1–2.4 (when using the area of the corresponding wing and the instantaneous velocity as reference area and reference velocity, respectively). This value is approximately twice as large as the steady-state value measured on a dragonfly wing at *Re*=730–1890 [steady-state aerodynamic forces on the fore- and hindwings of the dragonfly *Sympetrum sanguineum* were measured in a wind tunnel by Wakeling and Ellington(1997a); the maximum resultant force coefficient, obtained at an angle of attack of ∼60°, was approximately 1.3].

There are two possible reasons for the large vertical force coefficients of the flapping wings: one is the unsteady flow effect; the other is the effect of interaction between the fore- and hindwings (in the steady-state wind-tunnel test, interaction between fore- and hindwings was not considered).

### The effect of interaction between the fore- and hindwings

In order to investigate the interference effect between the fore- and hindwings, we computed the flow around a single forewing (and also a single hindwing) performing the same flapping motion as above. Fig. 8A,B gives vertical force(*C*_{L,sf}) and thrust (*C*_{T,sf})coefficients of the single forewing, compared with *C*_{L,f}and *C*_{T,f}, respectively. The differences between *C*_{L,sf} and *C*_{L,f} and between *C*_{T,sf} and *C*_{T,f} show the interaction effect. A similar comparison for the hindwing is given in Fig. 8C,D. For both the fore-and hindwings, the vertical force coefficient on a single wing (i.e. without interaction) is a little larger than that with interaction. For the forewing,the interaction effect reduces the mean vertical force coefficient by 14% of that of the single wing; for the hindwing, the reduction is 16% of that of the single wing. The interaction effect is not very large and is detrimental to the vertical force generation.

### The unsteady flow effect

The above results show that the interaction effect between the fore- and hindwings is small and, moreover, is detrimental to the vertical force generation. Therefore, the large vertical force coefficients produced by the wings must be due to the unsteady flow effect. Here, the flow information is used to explain the unsteady aerodynamic force.

First, the case of the single wing is considered. Fig. 9 shows the iso-vorticity surface plots at various times during one cycle. In order to correlate force and flow information, we express time during a stroke cycle as a non-dimensional parameter, *t̂*, such that *t̂*=0 at the start of the cycle and *t̂*=1 at the end of the cycle. After the downstroke of the hindwing has just started(*t̂*=0.125; Fig. 9A), a starting vortex is generated near the trailing edge of the wing, and a leading edge vortex (LEV)is generated at the leading edge of the wing; the LEV and the starting vortex are connected by the tip vortices, forming a vortex ring. Through the downstroke (Fig. 9B,C), the vortex ring grows in size and moves downward. At stroke reversal (between *t̂*≈0.36 and *t̂*≈0.65), the wing rotates and the LEV is shed. During the upstroke, the wing almost does not produce any vorticity. The vortex ring produced during the downstroke is left below the stroke plane (Fig. 9D–F)and will convect downwards due to its self-induced velocity. The vortex ring contains a downward jet (see below). We thus see that, in each cycle, a new vortex ring carrying downward momentum is produced, resulting in an upward force. This qualitatively explains the unsteady vertical force production. Fig. 10 shows the velocity vectors projected in a vertical plane that is parallel to and 0.6*R*from the plane of symmetry of the insect. The downward jet is clearly seen.

Fig. 11 shows the iso-vorticity surface plots for the fore- and hindwings (in the first half of the cycle the hindwing is in its downstroke; in the second half of the cycle the forewing is in its downstroke). Similar to the case of the single wing,just after the start of the first half of the cycle, a new vortex ring is produced by the hindwing (Fig. 11A); this vortex ring grows in size and convects downwards(Fig. 11A–C). Similarly,just after the start of the second half of the cycle, a new vortex ring is produced by the forewing (Fig. 11D), which also grows in size and convects downwards as time increases. Fig. 12 gives the corresponding velocity vector plots. The qualitative explanation of the large unsteady forces on the fore- and hindwings is similar to that for the single wing.

On the basis of the above analysis of the aerodynamic force mechanism, we give a preliminary explanation for why the forewing–hindwing interaction is not strong and is detrimental. The new vortex ring, which is responsible for the large aerodynamic force on a wing, is generated by the rapid unsteady motion of the wing at a large angle of attack. As a result, the effect of the wake of the other wing is relatively small. Moreover, the wake of the other wing produces downwash velocity, resulting in the detrimental effects.

### Power requirements

As shown above, the computed vertical force is enough to support the insect weight and the horizontal force is approximately zero; i.e. the force balance conditions of hovering are satisfied. Here, we calculate the mechanical power output of the dragonfly. The mechanical power includes the aerodynamic power(work done against the aerodynamic torques) and the inertial power (work done against the torques due to accelerating the wing-mass).

*C*

_{Q,a,t}and

*C*

_{Q,a,r}, respectively) are defined as:

*Q*

_{a,t}and

*Q*

_{a,r}are the aerodynamic torques around the axis of azimuthal rotation (

*z*′axis) and the axis of pitching rotation, respectively.

*C*

_{Q,a,t}and

*C*

_{Q,a,r}are shown in Fig. 13A,B. As can be seen,

*C*

_{Q,a,t}is much larger than

*C*

_{Q,a,r}.

*C*

_{Q,i,t}) is defined as:

*C*

_{Q,i,t}is shown in Fig. 13C. The inertial torque for rotation cannot be calculated since the moment of inertial of wing-mass with respect to the axis of flip rotation is not available. Because most of the wing-mass is located near the axis of flip rotation, it is expected that the inertial torque for rotation is much smaller than that for translation. That is, both the aerodynamic and inertial torques for rotation might be much smaller than those for translation. In the present study, the aerodynamic and inertial torques for rotation are neglected in the power calculation.

*C*

_{p}), i.e. power non-dimensionalized by 0.5ρ

*U*

^{3}(

*S*

_{f}+

*S*

_{h}),is:

*C*

_{p}of the fore- and hindwings is shown in Fig. 14. In the figure,contributions to

*C*

_{p}by the aerodynamic and inertial torques (represented by

*C*

_{p,a}and

*C*

_{p,i},respectively) are also shown. For the forewing(Fig. 14A), the time course of

*C*

_{p}is similar to that of

*C*

_{p,a}in the downstroke and to that of

*C*

_{p,i}in the upstroke; i.e. the aerodynamic power dominates over the downstroke and the inertial power dominates over the upstroke. This is also true for the hindwing(Fig. 14B).

Integrating *C*_{p} over the part of a wingbeat cycle where it is positive gives the coefficient of positive work(

*C*

_{p}over the part of the cycle where it is negative gives the coefficient of `negative' work (

*C*

^{-}

_{W}) for`braking' the wing in this part of the cycle.

*C*

^{+}

_{W}and

*C*

^{-}

_{W}for the forewing are 8.33 and –2.16, respectively. For the hindwing, they are 8.93 and –1.14, respectively.

*P*

^{*}) is defined as the mean mechanical power over a flapping cycle divided by the mass of the insect,and it can be written as follows (Sun and Tang, 2002):

*C*

_{W,f}and

*C*

_{W,h}are the coefficients of work per cycle for the fore- and hindwings, respectively. When calculating

*C*

_{W,f}or

*C*

_{W,h}, one needs to consider how the negative work fits into the power budget. There are three possibilities (Weis-Fogh,1972; Ellington,1984). One is that the negative power is simply dissipated as heat and sound by some form of an end stop, then it can be ignored in the power budget. The second is that in the period of negative work, the excess energy can be stored by an elastic element, and this energy can then be released when the wing does positive work. The third is that the flight muscles do negative work (i.e. they are stretched while developing tension, instead of contracting as in `positive' work) but the negative work uses much less metabolic energy than an equivalent amount of positive work and, again, the negative power can be ignored in the power budget. That is, out of these three possibilities, two ways of computing

*C*

_{W,f}or

*C*

_{W,h}can be taken. One is neglecting the negative work, i.e.:

*P*

^{*}is 37 W kg

^{-1}(when equations 22 and 23 are used,

*P*

^{*}is 30 W kg

^{-1}).

## Discussion

### Comparison with previous two-dimensional results

Wang (2000) and Lan and Sun(2001c) have presented two-dimensional (2-D) computations based on wing kinematics similar to those used in this study. Wang(2000) investigated a single airfoil; Lan and Sun (2001c)investigated both a single airfoil and fore and aft airfoils. It is of interest to make comparisons between the present three-dimensional (3-D) and the previous 2-D results.

The

*u*

_{t}is used as reference velocity and the

*u*

_{t}is used as reference velocity, the

^{2}=1.97]; approximately the same

*C*

_{L}of the forewing or the hindwing is nearly identical to that of the airfoil (compare Fig. 6A with fig. 3 of Wang, 2000).

Lan and Sun's results for the fore and aft airfoils showed that the interaction effect decreased the vertical forces on the airfoils by approximately 22% compared with that of the single airfoil(Lan and Sun, 2001c). For the fore- and hindwings in the present study, the reduction is approximately 15%,showing that 3-D forewing–hindwing interaction is weaker than in the 2-D case.

### Aerodynamic force mechanism and forewing–hindwing interaction

Recent studies (e.g. Ellington et al.,1996; Dickinson et al.,1999; Wu and Sun,2004) have shown that the large unsteady aerodynamic forces on flapping model insect wings are mainly due to the attachment of an LEV or the delayed stall mechanism. This is also true for the fore- and hindwings in the present study. The LEV dose not shed before the end of the downstroke of the fore- or hindwing (Fig. 11). If the LEV sheds shortly after the start of the downstroke, the LEV would be very close to the starting vortex, and a vortex ring that carries a large downward momentum (i.e. the large aerodynamic forces) could not be produced. Generation of a vortex ring carrying large downward momentum is equivalent to the delayed stall mechanism.

Data presented in Fig. 8show that the forewing–hindwing interaction is not very strong and is detrimental. In obtaining these data, the wing kinematics observed for a dragonfly in hovering flight (e.g. 180° phase difference between the forewing and the hindwing; no incoming free-stream) have been used. Although some preliminary explanation has been given for this result, we cannot currently distinguish whether or not this result will exist when the phasing,the incoming flow condition, etc., are varied. Analysis based on flow simulations in which the wing kinematics and the flight velocity are systematically varied is needed.

### Power requirements compared with quasi-steady results and with Drosophila results

Wakeling and Ellington(1997b,c)computed the power requirements for the dragonfly *Sympetrum sanguineum*. In most cases they investigated, the dragonfly was in accelerating and/or climbing flight. Only one case is close to hovering(flight SSan 5.2); in this case, the flight speed is rather low (advance ration is approximately 0.1) and the resultant aerodynamic force is close to the insect weight (see fig. 7Dof Wakeling and Ellington,1997b; fig. 5 of Wakeling and Ellington,1997c). Their computed body-mass-specific aerodynamic power is 17.1 W kg^{-1} (see table 3 of Wakeling and Ellington, 1997c;note that we have converted the muscle-specific power given in the table to the body-mass-specific power), only approximately half the value calculated in the present study. Lehmann and Dickinson(1997) and Sun and Tang(2002), based on experimental and CFD studies, respectively, showed that for fruit flies, calculation by quasi-steady analysis might under-estimate the aerodynamic power by 50%. A similar result is seen for the hovering dragonflies.

*P*

^{*}for the dragonfly in the present study (37 W kg

^{-1}) is not very different from that computed for a fruit fly (30 W kg

^{-1}; Sun and Tang, 2002), even though their sizes are greatly different (the wing length of the fruit fly is 0.3 cm and that of the dragonfly is 4.7 cm). For the fruit fly, the mechanical power is mainly contributed by aerodynamic power(Sun and Tang, 2002). It is approximately the case with the dragonfly in the present study (see Fig. 14). From equation 15 of Sun and Tang(2002), the aerodynamic torque of a wing can be written as:

*d̄*is the mean drag of the wing;

*r̂*

_{d}is the radius of the first moment of the drag normalized by

*R*. When the majority of the power is due to aerodynamic torque,

*P*

^{*}can be approximated as:

*r̂*

_{d}for the two insects is not very different. Then,

*P*

^{*}depends mainly on

*n*Φ

*R*(half the mean tip speed). The dragonfly's

*R*is approximately 16 times that of the fruit fly; but its

*n*Φ (36 H×69°) is approximately 1/14 of that of the fruit fly (240 H×150°). This explains why

*P*

^{*}of the dragonfly is not very different from that of the fruit fly.

**List of symbols**

- c
mean chord length of forewing

*C*_{d,f}drag coefficient of forewing

*C*_{d,h}drag coefficient of hindwing

*C*_{l,f}lift coefficient of forewing

*C*_{l,h}lift coefficient of hindwing

*C*_{L}total vertical force coefficient

- \(C_{\mathrm{L}}\)
total mean vertical force coefficient

*C*_{L,f}vertical force coefficient of forewing

*C*_{L,h}vertical force coefficient of hindwing

*C*_{L,sf}vertical force coefficient of single forewing

*C*_{L,sh}vertical force coefficient of single hindwing

*C*_{p}coefficient of power

*C*_{p,a}coefficient of aerodynamic power

*C*_{p,i}coefficient of inertial power

*C*_{Q,a,r}coefficient of aerodynamic torque for rotation

*C*_{Q,a,t}coefficient of aerodynamic torque for translation

*C*_{Q,i,t}coefficient of inertial torque for translation

*C*_{T}total thrust coefficient

- \(C_{\mathrm{T}}\)
total mean thrust coefficient

*C*_{T,f}thrust coefficient of forewing

*C*_{T,h}thrust coefficient of hindwing

*C*_{T,sf}thrust coefficient of single forewing

*C*_{T,sh}thrust coefficient of single hindwing

*C*_{W,f}coefficient of work per cycle of forewing

*C*_{W,h}coefficient of work per cycle of hindwing

- \(C_{\mathrm{W}}^{+}\)
coefficient of positive work

- \(C_{\mathrm{W}}^{-}\)
coefficient of negative work

- d̄
mean drag of a wing

*d*_{f}drag of forewing

*d*_{h}drag of hindwing

- I
moment of inertial of wing-mass

*l*_{f}lift of forewing

*l*_{h}lift of hindwing

- \(L\)
total mean vertical force

*L*_{f}vertical force of forewing

*L*_{h}vertical force of hindwing

- m
mass of the insect

- n
flapping frequency

- O,o,o′
origins of the two inertial frames of reference and the non-inertial frame of reference

- p
non-dimensional fluid pressure

*P*^{*}body-mass-specific power

*Q*_{a,t}aerodynamic torques around the axis of azimuthal rotation(

*z*′ axis)*Q*_{a,r}aerodynamic torques around the axis of pitching rotation

- r
radial position along wing length

*r*_{2}radius of the second moment of wing area of forewing

*r̂*_{d}radius of the first moment of wing drag

- R
wing length

- Re
Reynolds number

*S*_{f}area of one wing (forewing)

*S*_{h}area of one wing (hindwing)

- t
time

- t̂
non-dimensional time (

*t̂*=0 and 1 at the start and end of a cycle, respectively) *T*_{f}thrust of forewing

*T*_{h}thrust of hindwing

- u,v,w
non-dimensional velocity components in

*x,y,z*directions,respectively *u*_{t}translational velocity of a wing

- \(u_{\mathrm{t}}^{+}\)
non-dimensional translational velocity of a wing

*u*_{θ}azimuthal velocity

- U
reference velocity

- X,Y,Z
coordinates in inertial frame of reference (

*Z*in vertical direction) - x,y,z
coordinates in inertial frame of reference (

*z*perpendicular to stroke plane) - x′,y′,z′
coordinates in non-inertial frame of reference

- α
geometric angle of attack

- αd
midstroke geometric angle of attack of downstroke

- αu
midstroke geometric angle of attack of upstroke

- \({\dot{{\alpha}}}\)
angular velocity of flip rotation

- \({\dot{{\alpha}}}_{0}^{+}\)
non-dimensional angular velocity of flip rotation

- \({\dot{{\alpha}}}_{0}^{+}\)
a constant

- β
stroke plane angle

- γ
phase angle of the translation of a wing

- Δτr
duration of wing rotation or flip duration (non-dimensional)

- υ
kinematic viscosity

- ξ,η,ζ
transformed coordinates

- ρ
density of fluid

- τ
non-dimensional time

- τr
time when pitching rotation starts (non-dimensional)

- τc
period of one flapping cycle (non-dimensional)

- φ
azimuthal or positional angle

- Φ
stroke amplitude

- \({\dot{{\phi}}}\)
angular velocity of azimuthal rotation

- \({\dot{{\phi}}}^{+}\)
non-dimensional angular velocity of azimuthal rotation

- \({\ddot{{\phi}}}^{+}\)
non-dimensional angular acceleration of azimuthal rotation

## Acknowledgements

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

## References

**Alexander, D. E.**(

**Dickinson, M. H., Lehman, F. O. and Sane, S. P.**(

**Ellington, C. P.**(

**Ellington, C. P., van den Berg, C. and Willmott, A. P.**(

**Freymuth, P.**(

**Hilgenstock, A.**(

**Lan, S. L. and Sun, M.**(

**Lan, S. L. and Sun, M.**(

**Lan, S. L. and Sun, M.**(

**Lehmann, F.-O. and Dickinson, H. D.**(

*Drosophila melanogaster.*

**Meakin, R.**(

**Norberg, R. A.**(

**Norberg, R. A.**(

*Aeschna juncea*L., kinematics and aerodynamics. In

**Reavis, M. A. and Luttges, M. W.**(

**Rogers, S. E. and Kwak, D.**(

**Rogers, S. E., Kwak, D. and Kiris, C.**(

**Rogers, S. E. and Pulliam, T. H.**(

**Saharon, D. and Luttges, M.**(

**Saharon, D. and Luttges, M.**(

**Sane, S. P. and Dickinson, M. H.**(

**Somps, C. and Luttges, M.**(

**Sun, M. and Tang, J.**(

*Drosophila virilis.*

**Sun, M. and Wu, J. H.**(

**Wakeling, J. M. and Ellington, C. P.**(

**Wakeling, J. M. and Ellington, C. P.**(

**Wakeling, J. M. and Ellington, C. P.**(

**Wang, Z. J.**(

**Weis-Fogh, T.**(

*Drosophila.*

**Wu, J. H. and Sun, M.**(