A computational study of the aerodynamic forces and power requirements of dragonfly (Aeschna juncea) hovering

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.

The fore- and hindwings of the dragonfly are approximated by two flat plates. The thickness of the model wings is 1% of c (where cis 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₂=∫_Sfr²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 zaxis (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.

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, vand p with respect to y are set to zero).

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

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

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.

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.

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 (⁠

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.

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

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,fand 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 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.6Rfrom 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.

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

Integrating C_p over the part of a wingbeat cycle where it is positive gives the coefficient of positive work(⁠

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

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.

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.

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.

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₂

  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

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