SUMMARY
When an insect hovers, the centre of mass of its body oscillates around a point in the air and its body angle oscillates around a mean value, because of the periodically varying aerodynamic and inertial forces of the flapping wings. In the present paper, hover flight including body oscillations is simulated by coupling the equations of motion with the Navier–Stokes equations. The equations are solved numerically; periodical solutions representing the hover flight are obtained by the shooting method. Two model insects are considered, a dronefly and a hawkmoth; the former has relatively high wingbeat frequency (n) and small wing mass to body mass ratio, whilst the latter has relatively low wingbeat frequency and large wing mass to body mass ratio. The main results are as follows. (i) The body mainly has a horizontal oscillation; oscillation in the vertical direction is about 1/6 of that in the horizontal direction and oscillation in pitch angle is relatively small. (ii) For the hawkmoth, the peaktopeak values of the horizontal velocity, displacement and pitch angle are 0.11U (U is the mean velocity at the radius of gyration of the wing), 0.22c=4 mm (c is the mean chord length) and 4 deg., respectively. For the dronefly, the corresponding values are 0.02U, 0.05c=0.15 mm and 0.3 deg., much smaller than those of the hawkmoth. (iii) The horizontal motion of the body decreases the relative velocity of the wings by a small amount. As a result, a larger angle of attack of the wing, and hence a larger drag to lift ratio or larger aerodynamic power, is required for hovering, compared with the case of neglecting body oscillations. For the hawkmoth, the angle of attack is about 3.5 deg. larger and the specific power about 9% larger than that in the case of neglecting the body oscillations; for the dronefly, the corresponding values are 0.7 deg. and 2%. (iv) The horizontal oscillation of the body consists of two parts; one (due to wing aerodynamic force) is proportional to 1/cn^{2} and the other (due to wing inertial force) is proportional to wing mass to body mass ratio. For many insects, the values of 1/cn^{2} and wing mass to body mass ratio are much smaller than those of the hawkmoth, and the effects of body oscillation would be rather small; thus it is reasonable to neglect the body oscillations in studying their aerodynamics.
INTRODUCTION
When an insect hovers, the centre of mass of its body oscillates around a point in the air and its body angle (pitch angle) oscillates around a mean value, because of the periodically varying aerodynamic and inertial forces of the flapping wings. Similar oscillations exist in other flight modes (e.g. forward flight). In previous experimental and computational fluid dynamics (CFD) studies on insect aerodynamics [e.g. experimental studies (Ellington et al., 1996; Dickinson et al., 1999; Sane and Dickinson, 2001; Usherwood and Ellington, 2002a; Usherwood and Ellington, 2002b); CFD studies (Liu and Kawachi, 1998; Wang, 2000; Ramamurti and Sandberg, 2002; Sun and Tang, 2002a; Sun and Tang, 2002b)], the body dynamics was not included and the oscillations of the body centre of mass and body angle were neglected (the body was assumed to be fixed at a point or to move at a constant velocity). For insects with a relatively high wingbeat frequency and a small ratio of wing mass to body mass, the oscillations of the body would be very small and could be neglected. But for insects with a relatively low wingbeat frequency and a large ratio of wing mass to body mass, how large would be the effect of the body oscillations on the aerodynamic forces and moments? Moreover, in the study of flight dynamics (dynamic stability and control) of insects, the `equilibrium flight' must be given first. With body oscillations, the state variables at equilibrium flight are periodical functions of time. These periodical functions must be obtained before stability and control analysis can be conducted. As a result, it is of great interest to obtain the periodical solutions of the body motion at equilibrium flight.
In the present study, we couple the equations of motion with the fluid flow equations (the Navier–Stokes equations) to simulate the hovering (including body oscillations) of some model insects, obtaining the periodical solutions of the equilibrium flight. The results can tell us how body oscillations affect the aerodynamic forces and moments. Moreover, the results can provide the periodical solutions of equilibrium flight, based on which the flight stability and control problems can be studied.
Two model insects are considered. One is a model dronefly and the other a model hawkmoth. The model dronefly may represent insects with a relatively high wingbeat frequency and a small wing mass to body mass ratio; the model hawkmoth may represent insects with a relatively low wingbeat frequency and a large wing mass to body mass ratio (for the dronefly, the wingbeat frequency is about 160 Hz and the wing mass to body mass ratio is about 1%; for the hawkmoth, the corresponding values are 26 Hz and 6%). Another reason for choosing these insects is that their morphological and wing kinematic data at hovering are available.
MATERIALS AND METHODS
Equations of motion of an insect
Equations of motion for a model insect with a body and a number of flapping rigid wings have been derived by Sun and colleagues (Sun et al., 2007). Three frames of reference are used in the derivation (Fig. 1): frame (x_{E},y_{E},z_{E}) is an inertial frame (a frame fixed on the earth); x_{E} is in the horizontal direction and z_{E} in the vertical direction. Frame (x_{b},y_{b},z_{b}) is a frame fixed on the insect body with its origin at the centre of mass of the wingless body; x_{b}, y_{b} and z_{b} are the longitudinal, transverse and dorsoventral body axes, respectively. Frame (x_{w},y_{w},z_{w}) is a frame fixed on an insect wing with its origin at the root of the wing; x_{w} is in the spanwise direction, pointing to the wing tip, and z_{w} is in the chordwise direction, pointing to the trailing edge of the wing. With respect to the inertial frame, the centre of mass of the wingless body moves at velocity ν_{cg}, the body rotates at angular velocity ω_{bd} (angular velocity of the wing relative to the body is prescribed). The equations of motion can be written as (see Appendix A): (1) where _{b}ν_{cg} and _{b}ω_{bd} represent the x_{b}, y_{b} and z_{b} components of ν_{cg} andω _{bd}, respectively; and F represents the inertial, aerodynamic and gravitational forces and moments of the body and wings (see Appendix A).
Some terms in the righthand side of Eqn 1 are dependent on the orientation of the insect body, and the equations that govern the Euler angle rates of the body are needed. For longitudinal motion, the differential equation of the body pitch angle is: (2) where θ_{b} is the pitch angle of the body (angle between x_{b} and the horizontal) and q_{b} is the pitch rate, i.e. the y_{b} component ofω _{bd}.
The system of equations, Eqns 4 and 5, is integrated using a modified Euler method. The aerodynamic force and moment terms in Eqn 1 must be obtained by solving the fluid flow equations, i.e. the Navier–Stokes equations (see below for a description of how the equations of motion are coupled with the Navier–Stokes equations in the solution process).
Fluid dynamics equations
The aerodynamic forces and moments in Eqn 1 are determined by the Navier–Stokes equations: (3a) (3b) where u is the fluid velocity, p the pressure, ρ the density, ν the kinematic viscosity, = the gradient operator and =^{2} the Laplacian operator.
To obtain the aerodynamic forces and moments, in principle we need to compute the flows around the wings and the body. But near hovering, the aerodynamic forces and moments of the body are negligibly small compared with those of the wings because the velocity of the body is very small. Therefore, we only need to compute the flows around the wings. We further assume that the contralateral wings do not interact aerodynamically and neither do the body and the wings. Sun and Yu showed that the left and right wings had negligible interaction except during `clap and fling' motion (Sun and Yu, 2006); Aono and colleagues showed that the interaction between wing and body was negligibly small: the aerodynamic force of the case with the body–wing interaction differs by less than 2% from that without body–wing interaction (Aono et al., 2008). Therefore, the above assumptions are reasonable. Thus in the present flow computation, the body is neglected and the flows around the left and right wings are computed separately.
The computational methods used to solve the Navier–Stokes equations are the same as those described by Sun and Tang (Sun and Tang, 2002a). The algorithm was based on the method of artificial compressibility. The computational grid has dimensions 93×109×78 in the normal direction, around the wing section and in the spanwise direction, respectively; the normal grid spacing at the wall was 0.0015 (portions of the grids are shown in Fig. 2). The outer boundary was set at 20 chord lengths from the wing. The nondimensional time step was 0.02 (nondimensionalized by c/U where c is the mean chord of the wing of the insect and U is the mean flapping speed of the wing; see below for the definition of U). A detailed study of the numerical variables such as grid size, domain size, time step, etc., was conducted and it was shown that the above values for the numerical variables were appropriate for the calculations. Boundary conditions are as follows. For the farfield boundary condition, at the inflow boundary, the velocity components are specified as freestream conditions (determined by flight speed), whereas pressure is extrapolated from the interior; at the outflow boundary, pressure is set equal to the freestream static pressure and the velocity is extrapolated from the interior. On the wing surfaces, impermeable wall and nonslip conditions are applied and the pressure is obtained through the normal component of the momentum equation written in the moving grid system.
The integration method of coupling the equations of motion with the Navier–Stokes equations
The equations of motion (Eqns 1 and 2) and the Navier–Stokes equations (Eqn 3) must be coupled in the solution process because the aerodynamic forces and moments in the equations of motion must be obtained from the solution of the Navier–Stokes equations and the boundary condition of the Navier–Stokes equations must be obtained from the motion of the insect determined by the equations of motion. Let u_{b} and w_{b} denote the x_{b} and z_{b} components of , respectively, and q_{b} denote the y_{b} component of (note that in hovering, flight is symmetrical about the longitudinal plane and other components of and ω_{bd} are zero). The integration process is as follows. Suppose that at time t_{n} the motion of the insect body (u_{b}, w_{b}, q_{b} andθ _{b}) are known (the wing motion with respect to the body is prescribed), the boundary condition of the Navier–Stokes equations would be known and the flow equations can be solved to provide the aerodynamic forces and moments at this time station; then Eqns 1 and 2 are marched to the next time station t_{n}_{+1}, using the Euler predictorcorrector method; the process is repeated in the next step, and so on (see below for how the calculation starts). The Euler predictorcorrector method is described below.
Let x represent the vector of unknown variables of the body motion, or the state variables [u_{b} w_{b} q_{b}]^{T} and f(x,t) represent the vector function on the right side of the system of Eqns 1 and 2. The equations of motion can be written as: (4) The Euler predictorcorrector method has two steps (e.g. Gerald, 1978). First step, estimating (predicting) the value of x at t_{n}_{+1} using the single Euler method (let denote the estimated value): (5) Second step, obtaining the corrected value of x at t_{n}_{+1}: (6) The time accuracy of the integration method is of second order.
The method of finding the periodic solution
As mentioned above, in the flight of an insect, because of the cyclic force and moment of the flapping wings, the velocity of the centre of mass and the pitch angle of the body oscillate around some mean values. For hovering flight, the mean velocity and mean pitch rate are zero, and the mean pitch angle is a constant. Therefore, we are looking for periodic solutions of Eqns 1 and 2 with zero mean horizontal and vertical velocities and pitch rate. If the periodic solutions were stable, a forward integration of the equations would converge to the solutions. But linear stability analysis showed that hovering and forward flight were unstable (Taylor and Thomas, 2003; Sun and Xiong, 2005; Sun et al., 2007). Therefore, in the present study, the shooting method is used to find the periodical solutions. This method is described by Ling (Ling, 1983) (see also Fehlberg, 1969). In the shooting method, the initial conditions are iterated, so that they coincide with the terminal conditions. The method is outlined here and a detailed description is given in Appendix B. The equations of motion (Eqn 4) and the periodical conditions are: (7a) (7b) where T is the period of the solution (in the present study, T is known and is the wingbeat period). A vector parameter, s, is introduced, and it is chosen as x(0). Eqn 7 can be written as: (8a) (8b) where g(s) is a vector function representing the dependence of the equations on s; x(0; s)=s. s is iterated using Newton's method until the boundary conditions, i.e. F(s)=0, are satisfied to an acceptable accuracy (see Appendix B for the detailed procedure).
How the calculation starts and the solution process
In the shooting method described in the preceding section and Appendix B, the equations of motion are repeatedly integrated from t=0 to t=T with different initial conditions (initial state variables of the insect's body, u_{b}, w_{b}, q_{b} and θ_{b}), until the state variables at t=0 are equal to those at t=T. During each of the iterations, the Navier–Stokes equations need to be solved only in one wingbeat period (from t=0 to t=T). Note that flow is unsteady and the solution of unsteady Navier–Stokes in the period from t=0 to t=T depends on the flow field before t=0; specifically, on the flow field at the time step before t=0 (i.e. at t=–Δt; Δt is the time step). Since the flow is periodical in time, the flow field at t=–Δt is the same as that at t=T–Δt. But the flow field at t=T–Δt is not known until the equations of motion and the Navier–Stokes equations are solved in the period from t=0 to t=T.
How, then, can we obtain the flow field at t=–Δt? Let us look at the flow problem in which the model insect starts to beat its wings at t=0 (the first downstroke is from t=0 to t=0.5T) without body motion (see next section for the description of the wing motion). At t=0, the flow velocity is zero everywhere. With the known flow field at t=0, the Navier–Stokes equations are solved to give flows and aerodynamic forces at later times. Fig. 3 shows the computed lift, drag and pitching moment coefficients of the wing (denoted by C_{L}, C_{D} and C_{M}, respectively; see below for the definitions of the coefficients). In the first halfstroke (t=0 to t=0.5T), the aerodynamic forces are relatively large (this is because there is no previous wake to produce the induced effects). After less than three cycles, the forces and moment become approximately periodic. Usherwood and Ellington (experiment with revolving wings), Ramamurti and Sandberg (CFD computation of 3dimensional flapping wings), and Wang (CFD computation of 2dimensional wing) also reported that after the start of the wingbeat, the forces became periodic in less than three cycles (Usherwood and Ellington, 2002a; Ramamurti and Sandberg, 2002; Wang et al., 2004). This shows that the periodic wake of a flapping wing is established in less than three cycles after the wing starts to beat.
If the body of the insect oscillates, the established periodic wake of the flapping wing would be different from that in the above case of a fixed body. But if the velocity of the body oscillation is much smaller than that of the tip speed of the flapping wing, the difference in the established wake between these two cases would be small. This is because when the body velocity is much smaller than the speed of wing flapping, the wing flapping would dominate the wake production. Thus, when the body oscillation is small, we could use the established wake of the case without body oscillation to approximate that of the case with body oscillation. Let us test this by the following two flow computations. In the first one (case α), the body oscillates in the horizontal direction sinusoidally for four cycles (t=0 to t=4T); the amplitude of the horizontal velocity is 0.06U (where U is the mean flapping velocity of the wing, see below for its definition). In the second computation (case β), from t=0 to t=3T–ϵT (where ϵ is a small number, here ϵ=0.1), the body is fixed; from t=3T–ϵT to t=3T, the body velocity changes smoothly from zero to that of case α at t=3T; from t=3T to t=4T, the body motion is the same as that of case α (see Fig. 4A,B). In both cases, the flapping motion of the wing is the same. Note that in the fourth cycle (t=3T to t=4T), both the wing flapping motion and the body motion are the same for the two cases. Fig. 4C compares the computed C_{L} values in the fourth cycle of the two cases; they are almost identical (the same is true for the computed C_{D} and C_{M}, not shown here; in the above computation, ϵ is 0.1; we tried with ϵ varying between 0.05 and 0.2 and no difference was observed in the results of the fourth period). This shows that replacing the wake generated in the first three cycles by the flapping wing of the case with a small body oscillation (case α) with that of the case without body oscillation (case β) has negligible effect on the later flow computation.
Now, let us go back to the question raised above: in the calculation of coupling the equations of motion with the Navier–Stokes equations, the calculation is carried out in one period, t=0 to t=T; how is the flow field at t<0 given? On the basis of the above discussion, we assume that when the body oscillations are small, one can use the flow field established by the insect without body oscillation to provide the required initial flow field data.
With the above assumption, the solution process for obtaining the periodic solutions is as follows:

Guess the initial values of the state variables of the body, i.e. u_{b}, w_{b}, q_{b} andθ _{b} at t=0. In the first iteration, we set them as: u_{b}=0, w_{b}=0, q_{b}=0 and θ_{b}=χ (where χ is the mean body angle, known from measured data). It should be noted that in the later iterations, generally, the initial values of u_{b}, w_{b} and q_{b} would not be zero and that of θ_{b} would not be χ. For convenience, we write the initial values as s_{1}, s_{2}, s_{3} and s_{4}.

Start the wingbeat at t=–3T, with the body fixed (u_{b}, w_{b}, q_{b}=0 andθ _{b}=χ), and at the same time start to compute the flow with the Navier–Stokes equations. From t=–0.1T to t=0, u_{b}, w_{b}, q_{b} and θ_{b} change smoothly from zero to their initial values, s_{1}, s_{2}, s_{3} and s_{4}, respectively.

At t=0, start the computation of coupling the equations of motion and the Navier–Stokes equations (the flow field at the time step before t=0, which is needed for solving the Navier–Stokes equations at t≥0, is provided by the flow solution in step 2), and obtain the solution in the period from t=0 to t=T.

Compare the values of u_{b}, w_{b}, q_{b} and θ_{b} at t=T with those at t=0, i.e. s_{1}, s_{2}, s_{3} and s_{4}. If they are respectively equal, the periodic solution is obtained and the calculation finished. If not, go to the next step.

This step includes six `substeps' to obtain the Jacobian matrix in Newton's method (see Appendix B) and the new initial values.

Change the initial value of u_{b} by a small quantity,Δ s_{1}, and obtain a set of substep initial values: s_{1}+ Δs_{1}, s_{2}, s_{3} and s_{4}. Using the set of initial values, conduct the computation in steps 2 and 3 and obtain F(s^{*,1}).

Change the initial values of w_{b} byΔ s_{2} and obtain a second set of substep initial values: s_{1}, s_{2}+Δ s_{2}, s_{3} and s_{4}. Similarly, obtain F(s^{*,2}).

Similarly, obtain F(s^{*,3}).

Similarly, obtain F(s^{*,4}).

Compute the Jacobian matrix.

Using the Jacobian matrix, compute the new initial values.

With the new initial values, go to step 2 for the next iteration.
The wings, the flapping motion and the morphological and flight data
The model wings used in the present study are flat plates with rounded leading edges; the thickness is 3% of the mean chord length. The platforms of the model wings of the dronefly and the hawkmoth are obtained from the measured data of Liu and Sun, and Ellington (Liu and Sun, 2008; Ellington, 1984a), respectively (Fig. 2). The flapping motion of a wing is assumed as follows. The motion consists of two parts: the translation (azimuthal rotation) and the rotation (flip rotation); the outofplane motion of the wing (deviation) is neglected. The time variation of the positional angle (ϕ) of the wing is approximated by the simple harmonic function: (9) where n is the wingbeat frequency, θ the mean stroke angle andΦ the stroke amplitude. The angle of attack of the wing (α) takes a constant value during the downstroke or upstroke translation (the constant value is denoted by α_{d} for the downstroke translation andα _{u} for the upstroke translation); around stroke reversal, the wing flips and α changes with time, also according to the simple harmonic function. The function representing the time variation of α during the supination at the mth cycle is: (10a) where Δt_{r} is the time duration of wing rotation during the stroke reversal and a is a constant: (10b) t_{1} is the time when the wing rotation starts: (10c) The expression of the pronation can be written in the same way. For the dronefly, Eqns 9 and 10 are good approximations to the measured data (Liu and Sun, 2008; Ellington, 1984b); whilst for the hawkmoth, ϕ(t) is a little different from the simple harmonic function (Willmott and Ellington, 1997a; Willmott and Ellington, 1997b). From Eqns 9 and 10, we see that to prescribe the flapping motion, the stroke amplitudeΦ , the wingbeat frequency n, the time duration of wing rotationΔ t_{r}, the angles of attack in the downstroke and upstroke translations α_{d} and α_{u}, and the mean positional angle ϕ must be given.
Morphological and wing motion parameters for the dronefly are taken from the data of Liu and Sun (Liu and Sun, 2008) and that for the hawkmoth from the data of Willmott and Ellington (Willmott and Ellington, 1997a; Willmott and Ellington, 1997b). It should be noted that moments and products of inertia of the wing are not given in these references and they are estimated using the data of Liu and Sun, Willmott and Ellington, and Ennos (Liu and Sun, 2008; Willmott and Ellington, 1997a; Willmott and Ellington, 1997b), according to the density distribution of a dronefly wing measured by Ennos (Ennos, 1989b).
The morphological data for the wings are listed in Table 1: m_{wg} is the mass of one wing; R is wing length; c is the mean chord length of the wing; S is the area of one wing; r_{2}, is the radius of the second moment of wing area; I_{x}_{,w}, I_{y}_{,w} and I_{z}_{,w} are moments of inertia of the wing about the x_{w}, y_{w} and z_{w} axes, respectively; and I_{xy}_{,w} is the product of inertia of the wing.
The morphological data for the bodies are listed in Table 2: m_{b} is the mass of the body; I_{y,b} is the moment of inertia of the body about the y_{b}axis; l_{b} is body length; l_{1} is the distance between the wingbase axis and the centre of mass of the body; and χ_{0} is the free body angle [the definitions of χ_{0} and l_{1} are given in Fig. 5, following Ellington (Ellington, 1984a)]. For longitudinal motion, the body morphological parameters needed in equations of motion are the body mass, moment of inertia, and the relative position of the wingbase axis and the centre of mass. With χ_{0} and l_{1} known, the relative position of the wingbase axis and the centre of mass of the body can be determined. In the papers by Ellington and by Liu and Sun, the centre of mass (or values of χ_{0} and l_{1}) was determined as follows: two pictures showing the lateral and dorsoventral views of the body were taken; the crosssection of the body was taken as an ellipse and a uniform density was assumed, and then the position of the centre of mass can be computed (see Ellington, 1984a).
The wing and body kinematic data are listed in Table 3: Φ is the stroke amplitude; n is the wingbeat frequency; Δt_{r} is the wing rotation time; β is the stroke plane angle; χ is the mean body angle; Re, is the the Reynolds number of a wing (the mean flapping velocity, U, is defined as U=2Φnr_{2} and the Reynolds number of a wing is defined as Re=cU/ν where ν is the kinematic viscosity of the air).
From Tables 1, 2, 3, we see that all needed morphological and wingkinematic parameters have been given, exceptα _{d}, α_{u} and . These three parameters are difficult to measure accurately and they will be obtained during the solution process, using the condition of equilibrium flight (see below).
RESULTS
Test of the equations of motion and the computation codes
The code used for the flow computation was tested by Wu and Sun (Wu and Sun, 2004) using experimental data of flapping and rotating model insect wings; and recently, it was further tested by Sun and Yu (Sun and Yu, 2006) using experimental data for a flinging wing pair.
Here, we test the correctness of the equations of motion and their solution code by the following approach. Suppose that the model insect flaps its wings in the vacuum space without gravity; then an analytical solution of the problem exists, i.e. the centre of the total mass of the insect remains at the same point in space. Fig. 6 shows the time histories of the horizontal velocity of the centre of the total mass, together with that of the centre of the body mass, for the model dronefly which beats its wings in a horizontal plane withα _{d}=α_{u}=90 deg., =0,Φ =107 deg. and n=164 Hz (the vertical velocity of the centre of the total mass and that of the centre of the body mass are both zero and are not shown here; the calculation was carried up to a time of 300 wingbeats, but only results in the first 30 cycles are shown here). It can be seen that the centre of the total mass remains motionless; the centre of the body mass has a horizontal oscillation, because the wings oscillate in a horizontal plane. Integrations using various time step values showed that when the nondimensional time step value (nondimensionalized by c/U) was 0.02 (about 400 time steps in one cycle), the results changed by less than 0.02% when the time step was further halved.
Obtaining the periodical solutions of hover flight
If all the wing motion parameters are known, the periodical solutions of the hover flight can be solved directly using the methods given in Materials and methods. But as aforementioned, the downstroke and upstroke angles of attack (α_{d} and α_{u}) and the mean stroke angle () of the wings are not known; they are to be determined by the hover flight conditions: the mean vertical and horizontal velocities of the centre of mass of the body are zero and the mean pitch angle of the body is a constant (or the mean vertical aerodynamic force of the wings balances the weight of the insect and the mean horizontal aerodynamic force and the aerodynamic pitch moment around the centre of mass of the body are zero).
We determine the values of α_{d}, α_{u} and and the periodical solutions of the hovering insects in the following two steps.
Step 1: α_{d}, α_{u} and in the case of no body oscillations
In our previous work on the hovering flight of insects (e.g. Sun and Tang, 2002b; Sun and Du, 2003), the body was assumed to be fixed in the air and α_{d},α _{u} and were determined by the force and moment equilibrium conditions: a set of values of α_{d},α _{u} and were guessed; the flow equations were solved and the mean vertical force (Z_{e}), mean horizontal force (X_{e}) and mean pitch moment (M_{e}) were calculated; if the magnitude of Z_{e} was not equal to the weight, or X_{e} or M_{e} was not equal to zero, να_{d},α _{u} and were adjusted; the calculations were repeated until the magnitude of Z_{e} was approximately equal to the weight, and X_{e} and M_{e} were approximately equal to zero. Experimental observations showed that for many insects, α_{d} andα _{u} were about 35 deg. and was not far from zero (Ellington, 1984b). These values were used as the starting values in the computations. Using starting values of α_{d}, α_{u} and that were not far from the real values can reduce the number of repeated computations; more importantly, there could be more than one set of α_{d}, α_{u} and satisfying the equilibrium conditions due to the nonlinearity in the aerodynamic force production, and by so doing, the calculation could generally converge to the realistic solution. Since body oscillation was not considered, values of α_{d},α _{u} and determined in this way must be different from those when the body is free to oscillate. But it is expected that the difference would not be very large. As will be seen below,α _{d}, α_{u} and values determined in this way are useful in determining the α_{d},α _{u} and values when the body is free to oscillate.
Using the above method, we found that without body oscillations, for the dronefly, when α_{d}=α_{u}=31.8 deg. and =2.8 deg., the mean vertical force of the wings approximately balanced the weight, and the horizontal force and pitch moment are approximately zero. For the hawkmoth, the corresponding values areα _{d}=38.3 deg., α_{u}=39.5 deg. and =11 deg.
Step 2: solutions with body oscillations
The α_{d}, α_{u} and values obtained above give forces and moments satisfying the equilibrium conditions in the case of a body held fixed. When the body is free to oscillate, these α_{d}, α_{u} and values would not give the correct forces and moments for hovering. But it is expected that these α_{d},α _{u} and values are not too far from the correct ones, and that when these values are used, the insect would have a nearhovering flight, i.e. would fly at a very small velocity. Based on the flight speed and direction, we can adjust the values of α_{d},α _{u} and until an approximate hovering flight, i.e. with mean horizontal and vertical velocities and mean pitch rate being approximately zero, is obtained. Note that the horizontal and vertical velocities of the centre of mass of the body (denoted as ẋ_{E,b} and ż_{E,b}, respectively) are related to u_{b}, w_{b} and θ_{b} as: (11a) (11b)
First, let us obtain the solutions of hovering in the dronefly. Letα _{d}=31.8 deg., α_{u}=31.8 deg. and =2.8 deg. (the values for hovering without body oscillations). Solving Eqns 7a and 7b by the procedure described in Materials and methods (the initial values used are u_{b}=0, w_{b}=0, q_{b}=0 and θ_{b}=χ=42 deg.) gives the periodical solutions shown in Fig. 7A. For the results in Fig. 7 and below, the horizontal and vertical velocities of the centre of mass of the body, instead of u_{b} and w_{b}, are shown and the results have been nondimensionalized (the velocities are nondimensionalized by U, and q_{b} is nondimensionalized by c/U); Δθ_{b} denotes θ_{b} minus its mean value; for a clear description of the results, the time in a cycle is expressed as a nondimensional parameter, t̂, such that t̂=0 at the start of the downstroke and t̂=1 at the end of the subsequent upstroke (the integration time step, nondimensionalized by c/U, is 0.02; it was found that further reducing the time step produced little change in the results).
From Fig. 7A, it can be seen that when the model dronefly uses the above α_{d},α _{u} and values (31.8 deg., 31.8 deg. and 2.8 deg., respectively), the mean of ż_{E,b} is not zero but about 0.03U (the means of ẋ_{E,b} and q_{b} are 0.007U and 3×10^{–6}, respectively, approximately zero). That is, the insect is in descending flight at a low speed. This means that the values of α_{d} andα _{u} are not large enough for hovering. The values ofα _{d} and α_{u} are changed to 32.8 deg., and the resulting solutions are shown in Fig. 7B. It is seen that the mean of ẋ_{E,b} becomes a small negative value (about –0.015U) and the insect is flying upwards at a very low speed, showing that the values of α_{d} andα _{u} are a little too large. When α_{d} andα _{u} values are changed to 32.5 deg. (the resulting solutions are shown in Fig. 7C), the mean of ż_{E,b} is 0.003U, approximately zero (the means of ẋ_{E,b} and q_{b} are –0.002U and 3.3×10^{–6}, respectively, also approximately zero); that is, the insect is hovering. Thus the periodic solutions of hovering in the model dronefly are obtained.
Next, we obtain the solutions of hovering in the model hawkmoth. Using the same approach, it is found that when α_{d}=42 deg.,α _{u}=42.3 deg. and =10.5 deg., periodical solutions with zero means in ẋ_{E,b}, ż_{E,b} and q_{b} are obtained. The solutions are shown in Fig. 8.
It should be noted that there could be more than one periodic solution for a given set of morphological and kinematic data, because of the nonlinearity of the equations of motion. On the basis of experimental observations (e.g. Ellington, 1984b; Liu and Sun, 2008), the body oscillations are very small. In our computations, the initial values of the body motion variables are set as zero, and hence are not far from the realistic solution. Therefore, it is expected that the iterative computation could converge to the realistic solution.
DISCUSSION
The body oscillations
From the results in Fig. 7C and Fig. 8 it can be seen that during hovering, the centre of mass of a body mainly has a horizontal oscillation (oscillation in the vertical direction is relatively small: the amplitude of ż_{E,b} is about 1/6 that of ẋ_{E,b}). There are two reasons for the oscillation in the horizontal direction being much larger than that in the vertical direction. One is that for an insect in normal hovering (hovering with a horizontal stroke plane), the inertial force of the wing is approximately in the horizontal plane. The other reason is that during a wingbeat cycle, the frequency of the vertical aerodynamic force is twice as large as that of the horizontal aerodynamic force and, furthermore, the amplitude of oscillation of the vertical force is smaller than that of the horizontal force; this is explained as follows. The lift and drag of a wing are perpendicular and parallel to the stroke plane, respectively (e.g. Liu and Sun, 2008). Thus the vertical force is mainly contributed by the lift of the wing and the horizontal force by the drag of the wing. As shown in Figs 9 and 10, the lift coefficient (approximately equal to the vertical force coefficient) oscillates around its mean twice in one wingbeat cycle, and the drag coefficient (approximately equal to the horizontal force coefficient) oscillates around its mean only once.
For the hawkmoth, the peaktopeak value of the horizontal velocity of the centre of mass of the body is about 0.11U (Fig. 8) and that of the horizontal displacement of the centre of mass of the body (computed using the data in Fig. 8) is about 0.22c≈4 mm; the peaktopeak value of the pitch angle is 4 deg. For the dronefly (Fig. 7), the corresponding values are 0.02U, 0.05c≈0.15 mm and 0.3 deg. The body oscillations in the hawkmoth are much larger than those of the dronefly. This is expected because the wingbeat frequency of the hawkmoth is lower and the wing mass (relative to the body mass) is larger than that of the dronefly.
As an example, let us explain why the horizontal displacement of the body centre of mass (nondimensionalized by c) of the hawkmoth (0.22) is about 4.5 times as large as that of the dronefly (0.05). The horizontal displacement (H) can be written in two parts: (12) where H_{a} is due to the aerodynamic force and H_{i} is due to the inertial force of the wings. H_{a} and H_{i} can be approximately expressed as: (13) (14) where γ is the phase angle between the inertial and aerodynamic forces of the wings. The aerodynamic force and the inertial force of the wings can be approximately expressed as F_{a}sin2πnt and F_{i}sin(2πnt+γ), respectively. For H_{a}, we have: (15) or (16) F_{a} is due to the drag of the wings, and m_{b}g, the body weight, is equal to the mean lift of the wing. For many insects, the lift to drag ratio of the wing does not vary greatly (e.g. Usherwood and Ellington, 2002a; Usherwood and Ellington, 2002b; Wu and Sun, 2004); therefore, we assume F_{a}/m_{b}g does not vary greatly and obtain: (17) For H_{i}, we have: (18) or (19) F_{i} is proportional to 2m_{wg}r_{m} where r_{m} is the radius of the first moment of wing mass (the distance between the wing root and the centre of mass of the wing) and Eqn 19 can be written as: (20) Assuming r_{m}/c does not vary greatly, we obtain: (21) where 2m_{wg}/m_{b} is the ratio of wing mass (mass of two wings) to body mass.
Based on the data in Tables 1, 2, 3, the ratio of wing mass to body mass and the value of 1/cn^{2} of the model hawkmoth are 4–6 times as large as that of the model dronefly. This approximately explains the above results.
Comparison with experimental observations
Hedrick and Daniel measured the body oscillations of a hawkmoth hovering in front of an artificial flower and found that the centre of mass was maintained in a sphere of 6.5 mm radius (Hedrick and Daniel, 2006). The order of magnitude of our computed displacement (about 4 mm) is in agreement with these data. Liu measured the body motion of a hovering dronefly in a flight chamber and found that the displacement of the centre of mass was less than 0.5 mm (Liu, 2008). Again, the order of magnitude of the computed displacement (0.15 mm) is in agreement with the data.
It should be noted that the experimental observations were made under free flight conditions which included disturbances and stabilization control of the insect, whilst the present computation was not done under the same conditions. Therefore, the above comparisons are qualitative, with only the order of magnitude of the data being compared.
The effects of the body oscillations
Effect on the wing angles of attack needed for hovering
Table 4 gives the values ofα _{d}, α_{u} and needed for hover flight when the body is free to oscillate, compared with those when the body is held fixed. For the hawkmoth (which has a relatively low wingbeat frequency and a large wing mass to body mass ratio), the angles of attack needed to hover when the body is free to oscillate is about 3.5 deg. larger than that when the body is fixed. For the dronefly (which has a relatively high wingbeat frequency and a small wing mass to body mass ratio), this quantity is only 0.7 deg.
Here an explanation is given for why a larger angle of attack is needed when the body is free to oscillate. As seen in Fig. 8 (and also in Fig. 7C), in most of the downstroke (t=0 to t=0.5), the body moves backwards (ẋ_{E,b} is negative), thus the relative velocity of the wing is smaller than the case without body oscillations; in most of the upstroke (t=0.5 to t=10), the body moves forward (ẋ_{E,b} is positive), thus the relative velocity of the wing is also smaller than the case without body oscillations. As a result, a larger angle of attack is needed to produce enough weightbalancing force.
Effect on the aerodynamic forces on the wing
Fig. 9 gives the time courses of the lift (C_{L}) and drag (C_{D}) coefficients of the wing in one cycle for the hovering dronefly, with and without body oscillations, the forces are nondimensionalized by 0.5πU^{2}S, where S is the area of one wing. Fig. 10 gives the results for the hawkmoth.
Let us look at the results of the hawkmoth first (Fig. 10). For the wing lift (Fig. 10A), in some parts of the cycle (t̂≈0 to 0.25; t̂≈0.5 to 0.7), C_{L} in the case with body oscillations is slightly larger than that without body oscillations, and in other parts of the cycle (t̂≈0.25 to 0.5; t̂≈0.7 to 1.0), C_{L} in the case with body oscillations is slightly smaller than that without body oscillations. Thus the mean vertical force (mainly contributed by C_{L}) is the same for the cases with and without body oscillations, as it should be for balancing the weight. For the wing drag (Fig. 10B), in the whole cycle, the magnitude of C_{D} in the case with body oscillations is larger (by about 10%) than that without body oscillations. This is because the angle of attack of the wing in the case with body oscillations is larger than that without body oscillations. The larger wing drag would result in a larger power requirement (see below).
For the dronefly (Fig. 9), the force coefficients have rather small differences between the cases with and without body oscillations. This is expected from the above results that the body oscillations of the dronefly are very small.
Effect on power requirements
Based on the aerodynamic forces and the kinematic and morphological data, the mechanical power required for the flight can be computed [for the computational method, see for example Fry et al., or Liu and Sun (Fry et al., 2005; Liu and Sun, 2008)]. Fig. 11 gives the instantaneous nondimensional aerodynamic power, C_{p,a}, and inertial power, C_{p,i}, of a wing in one cycle for the model dronefly (the power is nondimensionalized by 0.5ρU^{2}Sc), for the cases with and without body oscillations. Fig. 12 gives the results for the model hawkmoth.
From C_{p,a} and C_{p,i}, the massspecific power (the mean mechanical power over a wingbeat cycle divided by the mass of the insect) can be calculated. In the case of zero elastic energy storage, the specific power is denoted as , and in the case of 100% elastic energy storage, as . The calculated results are given in Table 5. For the model hawkmoth, when the body is free to oscillate, the specific power is about 9% larger than that in the case without body oscillations. For the model dronefly, the corresponding value is about 2%.
Experimental and computational studies neglecting body oscillations
In previous experimental and computational studies on insect aerodynamics (e.g. Ellington et al., 1996; Dickinson et al., 1999; Sane and Dickinson, 2001; Sun and Tang, 2002a; Sun and Tang, 2002b; Berman and Wang, 2007), body oscillations were neglected. The above results show that for the hawkmoth, which has a relatively low wingbeat frequency and a large ratio of wing mass to body mass, when the body oscillations are included, the aerodynamic forces and power requirements differ by about 10% from that in the case of neglecting the body oscillations. But for the dronefly, which has a relatively high wingbeat frequency and small ratio of wing mass to body mass, the differences are rather small (about 2%).
As shown above, the motion of the body centre of mass in the horizontal direction dominates the body oscillation; and the nondimensional horizontal displacement of the body centre of mass consisted of two parts, one (due to wing aerodynamic force) proportional to 1/cn^{2} and the other (due to wing inertial force) proportional to the ratio of wing mass to body mass. Ellington, Ennos, and Willmott and Ellington obtained morphological data and hovering kinematic data for many insects [more than 20 species from 6 orders (Ellington, 1984a; Ellington, 1984b; Ennos, 1989a; Willmott and Ellington, 1997a; Willmott and Ellington, 1997b)]. These data show that hawkmoths are among those insects which have the largest values of 1/cn^{2} and wing mass to body mass ratio (n is about 25 Hz and wing mass to body mass ratio is about 6%). For many insects, their n is several times larger and the wing mass to body mass ratio several times smaller than that of the hawkmoth. For those insects, the differences in aerodynamic forces and power requirements between the cases with and without body oscillations would be much less than 10%. As a result, it can be said that for many insects, when studying their hovering aerodynamic forces and power requirements, it is reasonable to neglect the body oscillations.
The present results are useful to flight stability and control analysis
Flight dynamic stability and control are of great importance in the study of the biomechanics of insect flight. Flight dynamic stability deals with the perturbed motion of a flying body about its equilibrium state following a disturbance; if the perturbed motion decreases with time and goes to zero, the flight is stable; otherwise, it is unstable or neutrally stable. Flight control has two aspects. One is stabilization control, which deals with applying active control to stabilize the flight that is unstable or weakly stable; the other is control to change from one equilibrium state to another. Both the stability analysis and the control analysis need to start from the equilibrium state. The present method has provided the equilibrium state for hovering flight; it can also be used for obtaining periodic solutions of other flight modes (e.g. forward flight). In addition to finding the periodic solutions of equilibrium state, the integration process of coupling the equations of motion with the Navier–Stokes equations can be used to study the stability and control problems. But it should be noted that in the stability and control problems, lateral and forward translation and rotations about three axes of the body would occur, and the 6 degrees of body dynamics need to be considered, and, furthermore, the aerodynamic interactions of wing–wing and wing–body could become important. The dynamic equations in the present study (Eqn 1) are of 6 degrees of freedom. However, the Navier–Stokes solution might need to be modified to include the wing–wing and wing–body interaction effects.
Recently, some studies on flight dynamic stability and control were conducted using a simplified approach (e.g. Taylor and Thomas, 2003; Sun and Xiong, 2005). Two basic assumptions were made in the approach. The first was the `rigid body' assumption; it assumed that the frequency of the perturbed motion of the body was much lower than the wingbeat frequency and effects of the flapping wings on the flight system could be represented by the wingbeatcycle average aerodynamic and inertial forces. The second assumption was that linearization of the equations of motion was applicable. The validity of those assumptions needs to be tested, and this can be done by methods which solve the equations of the full nonlinear systems.
Some discussion on the solution method and the effects of variations in wing kinematics
The validity of the solution method
A major assumption of the solution method is that the historical wake of the flapping wing before the period of solution (t=0 to t=T) can be approximated by that in the case without body oscillations.
Now that we have obtained the periodic solutions of the body motion, we can examine the validity of the assumption as follows. We make two flow computations. In the first one (case A), from t=–3T to t=T, the model hawkmoth oscillates according to the periodic solution functions, u_{b}(t), w_{b}(t), q_{b}(t) andθ _{b}(t), obtained above. In the second computation (case B), from t=–3T to t=–0.1T, the body is fixed [u_{b}(t)=w_{b}(t)=q_{b}(t)=0,θ _{b}(t)=χ]; from t=–0.1T to t=0, u_{b}(t), w_{b}(t), q_{b}(t) andθ _{b}(t) change smoothly to u_{b}(0), w_{b}(0), q_{b}(0) andθ _{b}(0) of case A, respectively; from t=0 to t=T, u_{b}(t), w_{b}(t), q_{b}(t) andθ _{b}(t) are the same as those of case A (Fig. 13). The flapping motion of the wing, started at t=–3T, is the same for the two cases. The validity of the assumption can be examined by comparing the computed aerodynamic forces in the cycle from t=0 to t=T of the two cases. Fig. 13E shows the comparison. The results of the two cases are almost identical (in Fig. 13, only the comparison of C_{L} is given; comparisons of C_{D} and C_{M} between the two cases are similar), showing that the assumption is reasonable.
A further test is done by another two computations (called case a and case b). In case a, the body motion is similar to that of the above case A, except that the amplitudes of the oscillations are doubled. In case b, from t=–3T to t=–0.1T, the body is fixed; from t=–0.1T to t=0, u_{b}(t), w_{b}(t), q_{b}(t) and θ_{b}(t) change smoothly to u_{b}(0), w_{b}(0), q_{b}(0) and θ_{b}(0) of case a, respectively; and from t=0 to t=T, u_{b}(t), w_{b}(t), q_{b}(t) andθ _{b}(t) are the same as those of case a. The computed aerodynamic forces in the cycle from t=0 to t=T are given in Fig. 14. Only a very small difference is seen between the two cases. This shows that even if the body oscillations are twice as large as those of the model hawkmoth, the assumption is still reasonable.
As discussed above, the hawkmoth has a relatively low wingbeat frequency and a large ratio of wing mass to body mass. If the assumption is valid for body oscillations twice as large as that of the hawkmoth, it is expected that it is valid for many insects; it should be noted that for insects with rather large body oscillations, e.g. many butterflies [their velocity of body oscillation is comparable to wing flapping velocity (see Dudley, 1990)], the assumption would not be valid.
Effects of variations in wing kinematics
Idealized representations of wing kinematics are used in the present study, e.g. the outofplane motion of wing (deviation) is not included and the positional angle of the wing is approximated by the simple harmonic function. Here we make an additional computation for the hawkmoth to examine the effects of the approximations of wing kinematics on the solutions, by using the realistic wing kinematics (suggested by a reviewer of this paper).
Let the above computation for the hawkmoth [moth M1 in Willmott and Ellington (Willmott and Ellington, 1997a; Willmott and Ellington, 1997b)] be known as case 1. There the positional angleϕ (t) is the simple harmonic function, the deviation angleθ (t) is zero and the wing rotation time during the stroke reversal is 0.3T. The results of case 1 have been given above (Fig. 8).
The additional computation is known as case 2, in which more realistic wing kinematics is to be used. For the hawkmoth [moth M1 in Willmott and Ellington (Willmott and Ellington, 1997b)], only the time course of α(t) is available in Willmott and Ellington's paper. The time courses of ϕ(t) andθ (t) for this moth are available in another paper of Ellington and colleagues (Liu et al., 1998); they are replotted in Fig. 15A (the curves were based on the first three Fourier series terms, which were fitted to the measurements by Willmott and Ellington). ϕ(t) andθ (t) in Fig. 15A are used in case 2. As for α(t), we still assume it is constant in the middownstroke (denoted as α_{d}) and in the midupstroke (denoted as α_{u}), and during stroke reversal it is modeled by Eqn 10. The rotation time, determined from Willmott and Ellington's data [see figure 10B in Willmott and Ellington (Willmott and Ellington, 1997a)], is 0.3T for the supination and 0.2T for the pronation. The values of α_{d}, α_{u} and are adjusted to ensure the force and moment balance. The time derivatives of ϕ, θ and α are shown in Fig. 15B.
The computed results of case 2, compared with those of case 1, are given in Fig. 16. It can be seen that the main motion, i.e. the relatively large horizontal oscillation (ẋ_{E,b}), is changed only a little by the wing kinematics change (Fig. 16A), whilst the relatively small motions, the vertical oscillation (ż_{E,b}) and the body pitch oscillation (q_{b}, Δθ_{b}) are changed considerably by the wing kinematics change (Fig. 16B,C,D).
Fig. 17Ai, Bi and Ci shows the time course of the horizontal and vertical forces and pitching moment acting on the insect (coefficients C_{H}, C_{V} and C_{M} denote the nondimensional horizontal and vertical forces and pitching moment, respectively). The forces and moment are contributed by the aerodynamic and inertial forces and moment of the flapping wings. Fig. 17Aii and Aiii, Bii and Biii and Cii and Ciii shows the contributions from the aerodynamic and inertial forces and moments of the flapping wings. It is seen that variation in wing kinematics affects the aerodynamic forces and moments, as well as the inertial forces and moments of the wings. From Fig. 17Ai, Bi and Ci, it can be seen that the horizontal force (C_{H}) is much less affected by the wing kinematics change than the vertical force (C_{V}) and pitching moment (C_{M}). This is because the amplitude of oscillation of C_{H} is much larger than those of C_{V} and C_{M}, hence C_{H} is less sensitive to the wing kinematics change than C_{V} and C_{M}. This explains why the main motion, i.e. the horizontal oscillation of the body, is much less affected by the wing kinematics change than the relatively small vertical oscillation and body pitch oscillation.
APPENDIX A
For any vector V, in frame (x_{E},y_{E},z_{E}) (see Fig. 1), we write: (A1) where Vx_{E}, Vy_{E} and Vz_{E} are the x_{E}, y_{E} and z_{E} components of V, respectively. Similarly, in frame (x_{b},y_{b},z_{b}) and frame (x_{w},y_{w},z_{w}): (A2) where Vx_{b}, Vy_{b} and Vz_{b} are the x_{b}, y_{b} and z_{b} components of V, respectively, and Vx_{w}, Vy_{w} and Vz_{w} are the x_{w}, y_{w} and z_{w} components of V, respectively. _{E}V, _{b}V and _{w}V are related by the following (Sun et al., 2007): (A3) Here , and are matrix of direction cosines. For example, is determined by ψ_{b}, θ_{b} and ϕ_{b}, the Euler angles relating frame (x_{b},y_{b},z_{b}) to frame (x_{E},y_{E},z_{E}); is determined by ψ, ν_{w} θ_{w} andϕ _{w}, the Euler angles relating frame (x_{w},y_{w},z_{w}) to frame (x_{b},y_{b},z_{b}) (see Sun et al., 2007). With respect to the inertial frame, the centre of mass of the wingless body moves at velocity v_{cg}, the body rotates at angular velocityω _{bd} and a wing rotates at angular velocityω _{wg} (since angular velocity of the wing relative to the body is prescribed, ω_{wg} is related to ω_{bd} asω _{wg}=ω_{bd}+ω_{wg,0}, whereω _{wg,0} is the angular velocity of the wing relative to the body). R_{cg} is the position vector of the centre of mass of the body; R_{h} is the vector from the body centre of mass to the root of a wing; R_{wg} is the vector from the root of a wing to the centre of mass of the wing (Fig. 1). Let F_{A} be the total aerodynamic force and M_{A} the total aerodynamic moment about the centre of mass of the body, m_{b} the mass of the insect body, m_{wg} the mass of a wing, I_{bd} the matrix of moments and products of inertia of the body, I_{wg} the matrix of moments and products of inertia of a wing, g the gravitational acceleration and t the time. The equations of motion [equations 6 and 7 in Sun et al. (Sun et al., 2007)] are as follows: (A4) (A5) where N is the number of wings (in the present study, N=2).
The unknown variables in the equations of motion (Eqns A4 and A5) are ν_{cg} andω _{bd} (note that as aforementioned, ω_{wg} depends on ω_{bd} when the flapping motion of the wings,ω _{wg,0}, is prescribed). Solving for dν_{cg}/dt and dω_{bd}/dt, Eqns A4 and A5 can be rewritten as: (A6) where: where: I is a unit matrix, and expressions of F_{1}, F_{2}, F_{3}, M_{1},..., M_{5} are: where symbol ∼ is defined as follows: let a=[a_{1} a_{2} a_{3}]^{T} be a vector, and the ã denote the following matrix (see Etkin and Reid, 1972):
APPENDIX B
To solve the periodic system, Eqn 7, a vector parameter, s=x(0), is introduced, and Eqn 7 can be written as: (B1a) (B1b) where g(s) is a vector function representing the dependence of the equations on s; x_{0}=x(0; s)=s and x_{T}=x(T; s). s is iterated until the boundary conditions, i.e. F(s)=0, are satisfied to an acceptable accuracy. The iteration uses Newton's method: (B2) where the superscript i represents the iteration number, and F′(s) is the Jacobian matrix of F(s): (B3) where s_{j} and F_{k} are the jth and kth elements of s and F, respectively. F′_{jk}(s^{i}) can be computed using deference formulas as follows. Let: (B4) where Δs^{i}_{j} is a small quantity. Let F_{k}(s^{*,}^{j}) be the kth element of F(s^{*,}^{j}). F′_{jk}(s^{i}) is computed by: (B5)
The solution process is as follows.

An initial value x0= u_{b} w_{b} q_{b} θ_{b}]_{0}^{T} is guessed, denoted as s or x(0)^{1}, and Eqn B1a is solved numerically (by the integration process described in Materials and methods) up to t=T, obtaining x(T)^{1}, and then obtaining:

Choose a small quantity , denoted asΔ x_{1}(0)^{1}, and use: as the initial value to solve Eqn B1a and obtain the solution at t=T; the solution is denoted as [x(T)]^{*,1}. Thus, we obtain F(s^{*,1})=s^{*,1}–[x(T)]^{*,1}. F(s^{*,2}), F(s^{*,3}) and F(s^{*,4}) are calculated similarly.

Compute F′_{jk}(s^{i}) using Eqn B5, i.e.: (B6)

The improved initial value, s^{2} or x(0)^{2}, is obtained as: (B7)

Using x(0)^{2} as the initial value, repeat steps 1–4. The above process is repeated until x(0)^{i} and x(T)^{i} are approximately equal.
LIST OF ABBREVIATIONS
 a
 constant
 c
 mean chord length
 C_{D}
 drag coefficient
 C_{H}
 horizontal force coefficient
 C_{L}
 lift coefficient
 C_{M}
 nondimensional pitching moment coefficient
 C_{p,a}
 aerodynamic power
 C_{p,i}
 inertial power
 C_{V}
 vertical force coefficient
 CFD
 computational fluid dynamics
 F
 force
 F_{a}
 aerodynamic force of wing
 F_{A}
 total aerodynamic force
 F_{i}
 inertial force of wing
 g
 gravitational acceleration
 H
 horizontal displacement of body centre of mass (H_{i}, inertial; H_{a}, aerodynamic)
 I_{bd}
 matrix of moments and products of inertia of the body
 I_{mg}
 matrix of moments and products of inertia of the wing
 I_{x}_{,w}
 moment of inertia of the wing about the x_{w} axis
 I_{y}_{,w}
 moment of inertia of the wing about the y_{w} axis
 I_{z}_{,w}
 moment of inertia of the wing about the z_{w} axis
 l_{1}
 distance between wingbase axis and body centre of mass
 l_{b}
 body length
 M_{A}
 total aerodynamic moment about centre of mass of body
 M_{e}
 mean pitch moment
 m_{b}
 body mass
 m_{wg}
 mass of one wing
 N
 number of wings
 n
 wingbeat frequency
 P
 pressure
 massspecific power, zero elastic energy
 massspecific power, 100% elastic energy storage
 q_{b}
 pitch rate
 R_{cg}
 position vector of body centre of mass
 R_{h}
 vector of body centre of mass to root of wing
 R_{wg}
 vector from root of wing to centre of mass of wing
 Re
 Reynolds number
 r_{2}
 radius of second moment of wing area
 r_{m}
 radius of first moment of wing mass
 S
 area of one wing
 s
 vector parameter x(0)
 s_{1}, s_{2}, s_{3}, s_{4}
 initial values of u_{b}, w_{b}, q_{b} and θ_{b}
 T
 wingbeat period
 t
 time
 t_{1}
 time when wing rotation starts
 t̂
 nondimensional time
 u
 fluid velocity
 U
 mean velocity at radius of gyration of wing
 V
 vector
 X_{e}
 mean horizontal force
 x
 vector of unknown variables of body motion
 (x_{b},y_{b},z_{b})
 frame fixed on longitudinal body axis
 (x_{E},y_{E},z_{E})
 inertial frame in horizontal direction
 ẋ _{E,b}
 horizontal velocity of body centre of mass
 (x_{w},y_{w},z_{w})
 frame fixed on wing spanwise
 Z_{e}
 mean vertical force
 z_{b}
 frame fixed on dorsoventral body axis
 ż_{E,b}
 vertical velocity of body centre of mass
 gradient operator
 Laplacian operator
 α
 wing angle of attack (α_{d}, downstroke;α _{u}, upstroke)
 β
 stroke plane angle
 γ
 phase angle between inertial and aerodynamic forces of wings
 Δt_{r}
 duration of wing rotation
 θb
 pitch angle of body
 ν
 kinematic viscosity
 νcg
 velocity of centre of mass of wingless body
 ρ
 density
 Φ
 stroke amplitude
 ϕ
 positional angle of the wing
 ϕ
 mean stroke angle
 χ
 mean body angle
 χ0
 free body angle
 ωbd
 angular velocity of body
 ωwg
 angular velocity of wing
FOOTNOTES

This research was supported by grants from the National Natural Science Foundation of China (10732030), the 111 Project (B07009) and the Specialized Research Fund for the Doctoral Program of Higher Education (SRFDP, 200800061013).