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₂ (r₂=∫_Sf r²dS_f/S_f, where r is radial distance and S_f is the area of the forewing); r₂=0.61R (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₂ is called the translational speed.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 1.

Sketches of the model wings, the flapping motion and the reference frames. FW and HW denote fore- and hindwings, respectively. OXYZ is an inertial frame, with the X and Y axes in the horizontal plane; oxyz is another inertial frame, with the x and y axes in the stroke plane; o′x′y′z′ is a frame fixed on the wing, with the x′ axis along the wing chord and y′ axis along the wing span. β, stroke plane angle;φ , positional angle; α, angle of attack; R, wing length.

The flapping motion of a wing is simplified as follows. The wing translates downward and upward along the stroke plane and rotates during stroke reversal (Fig. 1B). The translational velocity is denoted by u_t and is given by: 1 where the non-dimensional translational velocity 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₂, where Φ and n are the stroke amplitude and stroke frequency of the forewing, respectively. Denoting the azimuth-rotational velocity as , we have .

The angle of attack of the wing is denoted by α. It assumes a constant value in the middle portion of a half-stroke. The constant value is denoted by α_d for the downstroke and by α_u for the upstroke. Around the stroke reversal, α changes with time and the angular velocity () is given by: 2 whereτ _r≤τ≤τ_r+Δτ_r, and the non-dimensional form ; is a constant; τ_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, can be determined (around the next stroke reversal, the wing would rotate fromα =180°–α_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

The Navier–Stokes equations are numerically solved using moving overset grids. For flow past a body in arbitrary motion, the governing equations can be cast in an inertial frame of reference using a general time-dependent coordinate transformation to account for the motion of the body. The non-dimensionalized three-dimensional incompressible unsteady Navier–Stokes equations, written in the inertial coordinate system oxyz (Fig. 1C), are as follows: 3 4 5 6 where 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.

The equations are first transformed from the Cartesian coordinate system (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: 7 Using equation 7, the transformation metrics in the inertial coordinate system, (ξ_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.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

Some portions of the moving overset grids.

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

The lift of a wing is the component of the aerodynamic force on the wing that is perpendicular to the translational velocity of the wing (i.e. perpendicular to the stroke plane); the drag of a wing is the component that is parallel to the translational velocity. 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: 8 9 These two formulae also apply to the case of the hindwing. The coefficients of 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: 10 etc., where ρ is the fluid density, 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: 11 12

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.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 was 10 000–30 000 deg. s^-1. Here, is set as 20 000 deg. s^-1, giving and Δτ_r=3.36.

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³. Fig. 3A shows the numerical and analytical drag coefficients (C_d) vs non-dimensional time (τ_s) (C_d=drag/0.5ρU_s²πa²;τ _s=tU_s/2a). 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/2a (r is radial distance) with fixed azimuthal angle π/2. The numerical result agrees well with the analytical solution outside the boundary layer.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 3.

Comparison between numerical and analytical solutions of a starting sphere. (A) Drag coefficient (C_d) vs non-dimensional time (τ_s). (B) Azimuthal velocity (u_θ) vs non-dimensional radial distance (r/2a).

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 30a 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.4a from the sphere surface and the out-boundary of the background grid is 30a 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.0c from the wing. The background Cartesian grid has dimensions of 90×85×80 and the outer boundary is 20c 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.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 4.

Comparison of the calculated and measured lift and drag forces. The experimental data are reproduced from fig. 3C,D of Sane and Dickinson (2001). (A,B) The midstroke angle of attack is 50° and stroke amplitude is 60°; (C,D) the midstroke angle of attack is 50° and stroke amplitude is 180°.

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 2c from the wing surface and that of the background grid at about 40c 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.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 5.

Non-dimensional angular velocity of flip rotation () and azimuthal rotation () of (A) hindwing and (B) forewing; (C) time courses of total vertical force coefficient (C_L) and (D) total thrust coefficient (C_T) in one cycle.

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 () [and the mean thrust coefficient ()]: =1.35 and =0.02. The value of 1.35 gives a vertical force of 756 mg, approximately equal to the insect mass (754 mg). The computed mean thrust (11 mg) is close to zero. That is, the force balance condition is approximately satisfied. In the calculation, the stroke plane angle, the midstroke angles of attack for the downstroke and the upstroke have been set as β=52°, α_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.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 6.

(A) Time courses of vertical force coefficients of forewing (C_L,f) and hindwing (C_L,h) and (B) thrust coefficients of the forewing (C_T,f) and the hindwing (C_T,h) in one cycle.

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.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 7.

(A) Time courses of lift coefficients of forewing (C_l,f) and hindwing (C_l,h) and (B) drag coefficients of the forewing (C_d,f) and the hindwing (C_d,h) in one cycle.

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²/(πU/2)²=2.1. Similarly, for the forewing, the peak value would be 2.6×[(S_f+S_h)/S_f]×U²/(πU/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.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 8.

(A) Time courses of vertical force coefficients of forewing (C_L,f) and single forewing (C_L,sf); (B) thrust coefficients of the forewing (C_T,f) and single forewing (C_T,sf); (C) vertical force coefficients of the hindwing (C_L,h) and single hindwing (C_L,sh) and (D) thrust coefficients of the hindwing (C_T,h) and single hindwing (C_T,sh) in one cycle.

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.6R from the plane of symmetry of the insect. The downward jet is clearly seen.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 9.

(A–F) Iso-vorticity surface plots at various times in one cycle (single hindwing). Note that the X axis is along the body of the dragonfly and the XZ plane is the plane of symmetry of the insect. t̂, non-dimensional time. The magnitude of the non-dimensional vorticity is 1.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 10.

Velocity vectors in a vertical plane parallel to and 0.6R from the plane of symmetry at various times in one cycle (single hindwing). The horizontal arrow at the top left represents the reference velocity (U). t̂, non-dimensional time.

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.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 11.

Iso-vorticity surface plots at various times in one cycle (fore- and hindwings). Note that the X axis is along the body of the dragonfly, and the XZ plane is the plane of symmetry of the insect. t̂, non-dimensional time. The magnitude of the non-dimensional vorticity is 1.

• Download figure
• Open in new tab
• Download powerpoint

• Download figure
• Open in new tab
• Download powerpoint

Fig. 12.

Velocity vectors in a vertical plane parallel to and 0.6R from the plane of symmetry at various times in one cycle (fore- and hindwings). The horizontal arrow at the top left represents the reference velocity (U). t̂, non-dimensional time.

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).

As expressed in equation 20 of Sun and Tang (2002), the aerodynamic power consists of two parts, one due to the aerodynamic torque for translation and the other to the aerodynamic torque for rotation. The coefficients of these two torques (denoted by C_Q,a,t and C_Q,a,r, respectively) are defined as: 13 14 where 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.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 13.

Time courses of aerodynamic torque coefficients for translation (C_Q,a,t) and rotation (C_Q,a,r) of (A) forewing and (B) hindwing in one cycle; (C) time courses of inertial torque coefficient for translation (C_Q,i,t) in one cycle.

The inertial power also consists of two parts (see equation 35 of Sun and Tang, 2002): one due to the inertial torque for translation and the other to the inertial torque for rotation. The coefficient of inertial torque for translation (C_Q,i,t) is defined as: 15 where is the non-dimensional angular acceleration of wing translation. 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.

The power coefficient (C_p), i.e. power non-dimensionalized by 0.5ρU³(S_f+S_h), is: 16 where 17 18 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).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 14.

Time courses of power coefficients of forewing (A) and hindwing (B) in one cycle. C_p, power coefficient; C_p,a, coefficient of power due to aerodynamic force; C_p,i, coefficient of power due to inertial force.

Integrating C_p over the part of a wingbeat cycle where it is positive gives the coefficient of positive work () for translation. Integrating C_p over the part of the cycle where it is negative gives the coefficient of `negative' work (C<sup>-</sup><sub>W</sub>) 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.

The body-mass-specific power (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): 19 where 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.: 20 21 The other is assuming the negative work can be stored and released when the wing does positive work, i.e.: 22 23 Here, equations 20 and 21 are used, and the computed P^* is 37 W kg^-1 (when equations 22 and 23 are used, P^* is 30 W kg^-1).

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 value (single airfoil) computed by Wang (2000) is approximately 1.97 [in fig. 4 of Wang (2000), maximum of u_t is used as reference velocity and the value is approximately 0.8; if the mean of u_t is used as reference velocity, the value becomes 0.8×(0.5π)²=1.97]; approximately the same value (single airfoil) was obtained by Lan and Sun (2001c). In the present study, the values for the single forewing and single hindwing are 1.51 and 1.64, respectively, approximately 20% less than the 2-D value. This shows that the 3-D effect on is significant. The wing length-to-chord ratio is not small (approximately 5); one might expect a small 3-D effect. But for a flapping wing (especially in hover mode), the relative velocity varies along the wing span, from zero at the wing base to its maximum at the wing tip, which can increase the 3-D effect. Note that although is reduced by the 3-D effect significantly, the time course of 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. 8 show 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. 7D of 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.

It is of interest to note that the value of 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: 24 where 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: 25 where is the ratio of the mean drag to the mean vertical force of the wing. For the fruit fly, this ratio is around 1 (Sun and Tang, 2002). For the dragonfly in the present study, since a large part of the vertical force is contributed by the drag, this ratio is not very different from 1. We assume that 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

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

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

coefficient of positive work

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

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₂

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

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

angular velocity of flip rotation

non-dimensional angular velocity of flip rotation

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

angular velocity of azimuthal rotation

non-dimensional angular velocity of azimuthal rotation

non-dimensional angular acceleration of azimuthal rotation