## Abstract

Because of the periodically varying aerodynamic and inertial forces of the flapping wings, a hovering or constant-speed flying insect is a cyclically forcing system, and, generally, the flight is not in a fixed-point equilibrium, but in a cyclic-motion equilibrium. Current stability theory of insect flight is based on the averaged model and treats the flight as a fixed-point equilibrium. In the present study, we treated the flight as a cyclic-motion equilibrium and used the Floquet theory to analyse the longitudinal stability of insect flight. Two hovering model insects were considered—a dronefly and a hawkmoth. The former had relatively high wingbeat frequency and small wing-mass to body-mass ratio, and hence very small amplitude of body oscillation; while the latter had relatively low wingbeat frequency and large wing-mass to body-mass ratio, and hence relatively large amplitude of body oscillation. For comparison, analysis using the averaged-model theory (fixed-point stability analysis) was also made. Results of both the cyclic-motion stability analysis and the fixed-point stability analysis were tested by numerical simulation using complete equations of motion coupled with the Navier–Stokes equations. The Floquet theory (cyclic-motion stability analysis) agreed well with the simulation for both the model dronefly and the model hawkmoth; but the averaged-model theory gave good results only for the dronefly. Thus, for an insect with relatively large body oscillation at wingbeat frequency, cyclic-motion stability analysis is required, and for their control analysis, the existing well-developed control theories for systems of fixed-point equilibrium are no longer applicable and new methods that take the cyclic variation of the flight dynamics into account are needed.

## 1. Introduction

Flight dynamic stability (inherent stability) is of great importance in the study of the biomechanics of insect flight. It is the basis for studying flight control, because the inherent stability of a flying system represents the dynamic properties of the basic system, such as which degrees of freedom are unstable, how fast the instability develops, and so on.

There has been some work on the flight dynamic stability of insects and robotic insects [1–5]. Taylor & Thomas [1] studied forward flight dynamic stability in desert locusts, providing the first formal quantitative analysis of dynamic stability in a flying animal. In the study, they used an averaged model and linear theory, borrowed directly from the literature of aircraft flight dynamics. They obtained the aerodynamic derivatives by measuring the aerodynamic forces and moments of tethered locusts in a wind tunnel at various wind speeds and body attitudes. Sun & Xiong [3] used the same theoretical framework to study the flight stability of a bumble-bee at hovering. They obtained the aerodynamic derivatives using the method of computational fluid dynamics (CFD; the computational approach allows simulation of the inherent stability of a flapping motion in the absence of active control, which is difficult to achieve in experiments using real insects). Sun and colleagues further studied the hovering stability in several insects with various sizes and wingbeat frequencies (hoverfly, dronefly, cranefly and hawkmoth) using the same methods [6,7]. Schenato *et al*. [2] and Deng *et al*. [4], based on the averaged model, developed stabilization control algorithms for the flapping flight of robotic insects.

Because of the periodically varying aerodynamic and inertial forces of the flapping wings, a flying insect is a cyclically forcing system, and generally it is not in a fixed-point equilibrium. For example, the centre of mass of a hovering insect oscillates around a spatial point and its body-pitch angle oscillates around some mean value. The ‘equilibrium flight’ is represented by a periodical solution of the equations of motion [8]. Thus, when studying the flight stability of an insect, in principle, one should consider the stability of the periodical solution that represents the ‘equilibrium flight’. But in the averaged model used in the above studies [1–4,6,7], the aerodynamic forces and moments were time-averaged over the wingbeat cycle, the equilibrium flight was represented by a fixed-point equilibrium and stability with respect to this equilibrium point was considered. Taylor & Zbikowski [5] studied the stability of forward-flying locusts without assuming the equilibrium as a fixed point. They used instantaneous force and moment measurements from tethered locusts to parametrize the equations of motion and developed a nonlinear time-periodical model of the locusts. By numerically integrating the nonlinear equations in the model, they explored the stability problem of the locusts. Because forces and moments measured from real insects were used, some control responses were necessarily included in the model. As a result, the results they obtained were that for a controlled system; the model could not provide information on the inherent stability of an insect. Moreover, by numerical simulation, the study only gave specific solutions for the given initial data, while the general stability property could not be obtained.

In the present study, we analyse the stability of the periodic solution that represents the ‘equilibrium flight’, i.e. the cyclic-motion stability, using the Floquet theory (when a periodic solution is slightly perturbed, if the perturbation dies out as time increases, the periodic solution is stable, otherwise it is unstable; the equations governing the small perturbation are a system of linear ordinary differential equations with time-periodic coefficients, and the theory that has been developed to solve this type of equation is known as Floquet theory). The Floquet stability analysis can give the natural modes of motion of the perturbation, and hence the general stability property of the cyclic-motion equilibrium. It is expected that the results of the Floquet theory (cyclic-motion stability analysis) would be close to that of the averaged model (fixed-point stability analysis) for small insects, which generally have very small amplitude of body oscillation because their wingbeat frequency is relatively high and wing-mass to body-mass ratio being relatively small (analysis by Wu *et al*. showed that the amplitude of body oscillation consisted of two parts: one, owing to wing aerodynamic force, was proportional to 1/*cn*^{2} (*c* is the mean chord length of wing and *n* the wingbeat frequency) and the other, owing to wing inertial force, was proportional to wing-mass to body-mass ratio [8]), but for relatively large insects, which have relatively large amplitude of body oscillation, the Floquet theory could give better results.

We consider the dynamic stability of hover flight in two representative model insects—a model dronefly and a model hawkmoth; the former may represent insects with small amplitude of body oscillation and the latter may represent insects with relatively large amplitude of body oscillation (the peak-to-peak displacement of the body's centre of mass of a hovering dronefly is only about 1% of the body length, but that of a hovering hawkmoth is about 10% of the body length [8–10]). As a first step, we study the longitudinal flight stability. First, we obtain the periodic solution of the equilibrium flight by numerically solving the equations of motion coupled with the Navier–Stokes equations. Then, we study the stability of the periodic solution using the Floquet theory.

## 2. Material and methods

### 2.1. Equations of motion

Let (*x*_{b}, *y*_{b}, *z*_{b}) be a frame fixed on the insect body with its origin at the centre of mass of the wingless body, *x*_{b} and *z*_{b} be in the longitudinal symmetrical plane (pointing forward and downward, respectively) and *y*_{b} points to the right side of the insect (figure 1). The longitudinal equations of motion are (see electronic supplemental material for the derivation of the equations; note that the electronic supplementary material contains some important details of the study):
2.1a
2.1b
2.1c
2.1d

where *u* and *w* denote the *x*_{b} and *z*_{b} components of velocity of the insect body, respectively, and *q* denotes the *y*_{b} component of the angular velocity of the body; *m* the total mass of the insect (body mass plus wing mass); *θ* the angle between *x*_{b} and the horizontal (in the present study, *x*_{b} is so chosen that during hovering its mean direction is horizontal, and thus during hovering, the mean of *θ* is zero) and *I*_{y,bw} is the moment of inertia of the body and wings about *y*_{b} axis; *XA* and *ZA* are the *x*_{b} and *z*_{b} components of the total aerodynamic force, respectively, and *MA* is the *y*_{b} component of the total aerodynamic moment about the centre of mass of the body; *XI* and *ZI* are the *x*_{b} and *z*_{b} components of the inertial forces owing to flapping wings, respectively, and *MI* is the *y*_{b} component of the moment owing to the wing inertial forces (*XI* and *ZI* are called wing-inertial forces and *MI* is called the wing-inertial moment); *MG* is the *y*_{b} component of the moment produced by the weight of the wings; *u*, *w*, *q* and *θ*, which determine the motion of the flying system, are called state variables.

### 2.2. Navier–Stokes equations and their solution method

The Navier–Stokes equations that determine the aerodynamic force and moment in equation (2.1) are: 2.2aand 2.2b

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 when compared with those of the wings because the velocity of the body motion is very small, and 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 & Yu [11] showed that the left and right wings had negligible interaction except during ‘clap and fling’ motion; Yu and Sun [12] showed that interaction between wing and body was negligibly small: the aerodynamic force in the case with body–wing interaction is less than 2 per cent different from that without body–wing interaction. Therefore, the above assumptions are reasonable. Thus, in the present CFD model, the body is neglected and the flows around the left and right wings are computed separately.

The computational method used to solve the Navier–Stokes equations was the same as that described by Sun & Tang [13]. It was based on the method of artificial compressibility and was developed by Rogers *et al.* [14]. In the method, the time derivatives of the momentum equations were 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 was introduced into the equations and a pseudo-time derivative of pressure divided by an artificial compressibility constant was introduced into the continuity equation. The resulting system of equations was iterated in pseudo-time until the pseudo-time derivative of pressure approached zero, thus, the divergence of the velocity at the new time level approached zero. The derivatives of the viscous fluxes in the momentum equation were approximated using second-order central differences. For the derivatives of convective fluxes, upwind differencing based on the flux-difference splitting technique was used. A third-order upwind differencing was used at the interior points and a second-order upwind differencing was used at points next to boundaries. Details of this algorithm can be found in the paper by Rogers *et al.* [14].

The computational grids were generated using a Poisson solver and they were O–H type grids (the grids are further described in the electronic supplementary material and the results of studies on the convergence of solutions are also given there).

Boundary conditions were as follows. For the far-field boundary condition, at the inflow boundary, the velocity components were specified as free-stream conditions (determined by known wing motion), whereas pressure was extrapolated from the interior; at the outflow boundary, pressure was set equal to the free-stream static pressure and the velocity was extrapolated from the interior. On the wing surface, impermeable wall and non-slip conditions were applied and the pressure was obtained through the normal component of the momentum equation written in the moving grid system.

### 2.3. Wings, flapping motion and flight data

Using measured wing deformation data [15], Du & Sun [16] investigated the effect of wing deformation on aerodynamic forces in hovering hoverflies and showed that as a first approximation, the deformable wing could be modelled by a rigid flat-plate wing with its angle of attack being equal to the local angle of attack at the radius of the second moment of wing area. Thus, in the CFD model, we further assumed that wings were rigid flat plates. Hawkmoth wings have relatively large deformation during flapping flight [17,18] and the rigid-wing model could be less accurate for hawkmoths. But as the main purpose of the present work was to examine the difference between Floquet and averaged-model theories, and the rigid-wing model was used in both theories, we expected that using the rigid-wing model was reasonable. The platforms of the model wings of the dronefly and the hawkmoth were obtained from the measured data of Liu & Sun and Ellington, respectively [10,19], and the wing section was a flat plate of 3 per cent thickness with rounded leading and trailing edges.

On the basis of experimental data [10,20–22], the flapping motion of a wing was assumed as follows. The motion consisted of two parts: the translation (azimuthal rotation) and the rotation (flip rotation) (the out-of-plane motion of the wing (deviation) is neglected). The time variation of the positional angle (*ϕ*) of the wing was approximated by the simple harmonic function:
2.3

where *n* was the wingbeat frequency, the mean stroke angle and *Φ* the stroke amplitude. The angle of attack of the wing (*α*) took a constant value during the down- or upstroke translation (the constant value was denoted by *α*_{d} for the downstroke translation and *α*_{u} for the upstroke translation); around stroke reversal, the wing flipped and *α* changed with time according to the simple harmonic function. The function representing the time variation of *α* during the supination at *m*th cycle was:
2.4

where Δ*t*_{r} was the time duration of wing rotation during the stroke reversal and *a* a constant
2.5a*t*_{1} was the time when the wing rotation started
2.5bwhere *T* was the wingbeat cycle. The expression of the pronation could be written in the same way.

For the dronefly, equations (2.3)–(2.5) are good approximations of the measured data [10,20]. For the hawkmoth, *ϕ*(*t*) is a little different from the simple harmonic function [21,22]. The effect of this difference on body motion was investigated by Wu *et al*. and it was shown that the motion was changed only a little by the difference (see fig. 16A of Wu *et al*. [8]). From equations (2.3)–(2.5), we see that to prescribe the flapping motion, *Φ*, *n*, Δ*t*_{r}, *α*_{d}, *α*_{u} and must be given.

Morphological and wing motion parameters for the dronefly were taken from the data by Liu & Sun [10] and that for the hawkmoth from the data by Willmott & Ellington [21,22]. It should be noted that moments and products of inertia of wing were not given in these references and they were estimated using the data on density distribution of a dronefly wing [23].

Morphological data for the wings are listed in table 1: *m*_{wg}, mass of one wing; *R*, wing length; *c*, mean chord length of wing; *S*, area of one wing; *r*_{2}, radius of second moment of wing area; *I*_{x,w}, *I*_{y,w} and *I*_{z,w}, moments of inertia of wing about *x*_{w}, *y*_{w} and *z*_{w} axes, respectively; *I*_{xy,w}, product of inertia of wing. Morphological data for the bodies are listed in table 2: *m*_{bd}, mass of body; *I*_{y,b}, moment of inertia of body about *y*_{b} axis; *l*_{b}, body length; *l*_{1}, distance between the wing-base axis and the centre of mass of the body; *χ*_{0}, free body angle (see the study of Ellington [19] for the definitions of *χ*_{0} and *l*_{1}; with χ_{0} and *l*_{1} known, the relative position of the wing-base axis and the centre of mass of the body can be determined). *I*_{y,bw} required in equation (2.1) is determined by *I*_{y,b} and *I*_{x,w}, *I*_{y,w}, *I*_{z,w}, and *I*_{xy,w}.

Wing and body kinematic data are listed in table 3: *Φ*, *n*, Δ*t*_{r}, stroke plane angle *β*, mean body angle *χ*, Reynolds number of wing *Re* (*Re* = *cU*/ν, where ν is the kinematic viscosity of the air and *U* the mean flapping velocity, *U* = 2*Φ**nr*_{2}).

From tables 1 to 3, we see that all needed wing-kinematic parameters are given, except *α*_{d}, *α*_{u} and . These three parameters are determined in the process of obtaining the periodical solution of hovering. The reasons for not using experimental data for the three parameters are as follows. There are necessarily small errors in the measurement. When there are small errors in the wing-kinematic parameters, a periodic solution not representing a hover flight, but a flight with some low speed is obtained [8]; for example, if the measured angle of attack is a little larger than the real one, the computed periodic solution gives an upward flight of low speed. To overcome this problem, we do not use experimental data to prescribe α_{d}, *α*_{u} and , but determine them in the solution process, by requiring that they take values such that hovering flight is ensured. The reason for only three parameters being determined in this way is that there are three conditions to be met in hovering flight: zero mean horizontal velocity, zero mean vertical velocity and zero mean pitch rate. The reason for choosing *α*_{d}, *α*_{u} and to be determined in the solution process is that experimental data for *α*_{d} and *α*_{u} have relatively large error [10], and aerodynamic forces and moments are very sensitive to the variations in *α*_{d}, *α*_{u} and .

### 2.4. The periodic solution

The motion of a flying insect is governed by the equations of motion (equation (2.1)) coupled with the Navier–Stokes equations (equation (2.2)). At constant-speed flight or hovering, the wings flap with constant kinematic parameters and the inertial (*XI*, *ZI* and *MI*) and aerodynamic (*XA*, *ZA* and *MA*) forces and moments are periodic. Thus, the system represented by equations (2.1) and (2.2) is a periodically driving system. Does the system have a periodic solution, or a fixed-point or non-periodical but oscillatory solution [4]? Because the system consists of highly nonlinear ordinary (equation (2.1)) and partial (equation (2.2)) differential equations, mathematically proving the existence of periodic solution is very difficult. Hovering insects are observed to oscillate about a fixed point (e.g. a hawkmoth hovering in front of an artificial flower [9]; a dronefly hovering in a flight chamber [24]). Here, on the basis of experimental observations, we assume that a periodic solution of equation (2.1) (coupled with equation (2.2)), which represents the hover flight, exists. We obtain the periodic solution by numerically solving equation (2.1) coupled with equation (2.2).

A method of obtaining the periodic solution of hover flight has been developed by Wu *et al*. [8]. A brief description of the method and the solution process is given in the electronic supplementary material.

### 2.5. Equations of disturbance motion

Let *u*_{0}(*t*), *w*_{0}(*t*), *q*_{0}(*t*) and *θ*_{0}(*t*) be the periodic solution representing the hover flight. They satisfy equation (2.1):
2.6a
2.6b
2.6c
2.6dwhere *XA*_{0}(*t*) and *ZA*_{0}(*t*) are the periodic aerodynamic forces and *MA*_{0}(*t*) the periodic aerodynamic moment at hover flight; *XI*_{0}(*t*), *ZI*_{0}(*t*) and *MI*_{0}(*t*) are the corresponding wing-inertial forces and moment. The period of the above periodic functions is *T*.

When the equilibrium flight is disturbed, we express each state variable as the sum of the equilibrium-flight value and the disturbance value:
2.7a
2.7b
2.7c
2.7dwhere the symbol *δ* denotes a small quantity. The aerodynamic forces and moment will change from the equilibrium value when the insect is in disturbance motion, and so will the inertial forces and moment of the wings because some parts of the inertial forces and moment (e.g. wing inertial force owing to gyroscopic effect of the wing and body rotations) are dependent on the body motion variables. We also express the aerodynamic and inertial forces and moments as the sum of the equilibrium flight value and the disturbance value:
2.8a
2.8b
2.8c
2.9a
2.9b
2.9c

Substituting equations (2.7)–(2.9) into equation (2.1), using equation (2.6), and neglecting the second and higher order terms give 2.10a 2.10b 2.10c 2.10d

Taylor's theorem for analytical functions is used to express each of the disturbance values of the aerodynamic and wing-inertial forces and moments as an infinite series about the reference condition, retaining only the first-order terms to give an expansion of the form:
2.11a
2.11b
2.11c
2.12a
2.12b
2.12cwhere *∂XA*(*t*)/*∂u*, *∂ZA*(*t*)/*∂u*, etc., are partial derivatives of the aerodynamic and wing-inertial forces and moments at equilibrium flight; their definition is given as follows. Let us take the derivatives with respect to *u*, i.e. *∂XA*(*t*)/*∂u*, *∂ZA*(*t*)/*∂u* and *∂MA*(*t*)/*∂u*, as an example. At equilibrium flight, as described above, the state variables, *u*_{0}(*t*), *w*_{0}(*t*), *q*_{0}(*t*) and *θ*_{0}(*t*), and the aerodynamic forces and moment, *XA*_{0}(*t*), *ZA*_{0}(*t*) and *MA*_{0}(*t*), are periodic functions of time with *T* as period. Let a small constant value of *u*, denoted as Δ*u*, be added to the hovering motion variables, i.e. the motion of the insect becomes *u*_{0}(*t*) + Δ*u*, *w*_{0}(*t*), *q*_{0}(*t*) and *θ*_{0}(*t*), and let the aerodynamic forces and moment at this motion be

The partial derivatives with respect to *u*, i.e. *∂XA*(*t*)/*∂u*, *∂ZA*(*t*)/*∂u* and *∂MA*(*t*)/*∂u*, are defined, respectively, as
2.13a
2.13b
2.13c

The derivatives with respect to *w* or *q*, and the inertial force and moment derivatives are defined in the same way. In addition to giving the definition of the partial derivatives, equation (2.13) also shows how the derivatives can be computed approximately.

In the above definition, the motion variables, *u*_{0}(*t*)+Δ*u*, *w*_{0}(*t*), etc., are periodic functions of time with period *T* because Δ*u* is constant, and hence Δ*XA*(*t*), Δ*ZA*(*t*) and Δ*MA*(*t*), or *∂XA*(*t*)/*∂u*, *∂ZA*(*t*)/*∂u* and *∂MA*(*t*)/*∂u* are periodic functions of time with the same period. Similarly, the other aerodynamic derivatives and the wing-inertial force derivatives are periodic functions of time with period *T*.

The disturbance values of aerodynamic forces and moments, *δ**XA*(*t*), *δ**ZA*(*t*), etc., are a function of the state variables and their time derivatives, , … (where ‘.’ denotes time differentiation). But in equations (2.11) and (2.12), we have only included terms owing to the variation of the state variables (*δ**u*, *δ**w* and *δ**q*), and terms due to the time derivatives of state variables (, etc.) have been neglected. That is, we have assumed that unsteady aerodynamic effects due to the disturbance motion of the insect body are much smaller than that due to the flapping motion of the wings. The validity of the assumption will be tested below by numerically solving the full equations of motion coupled with the Navier–Stokes equations.

Substituting equation (2.11) and (2.12) into equation (2.10) gives 2.14a

where **A**(*t*) is called as system matrix and

2.14b

Let us point out that in equation (2.14) and henceforth, all the variables are non-dimensionalized, using *c*, *U* and *T* (*T* = 1/*n*) as the reference length, velocity and time, respectively: time *t* is non-dimensionalized by *T*; a velocity by *U*; an angular velocity by 1/*T*; a force by 0.5*ρ*U^{2}(2*S*) (*S* is area of one wing); a moment by 0.5*ρ*U^{2}(2*S*)*c*; *m* by 0.5*ρ**U*(2*S*)*T*; a moment of inertia by 0.5*ρ**U*^{2}(2*S*)*cT ^{2}*;

*g*by

*U*/

*T*(

*ρ*is 1.25 kg m

^{−3}and

*g*is 9.81 m s

^{−2}).

Each element of **A**(*t*) is a periodic function of time with non-dimensional period 1, and equation (2.14) is a system of linear ordinary differential equations with periodic coefficients.

### 2.6. Floquet stability analysis of the periodic solution

Equation (2.14) can be solved to yield insights into the stability of the periodic solution that represents hover flight. If the perturbations, *δ**u*(*t*), *δ**w*(*t*), *δ**q*(*t*) and *δ**θ*(*t*), die out as time increases, the periodic solution or the hover flight is stable, otherwise it is unstable. The general theory for such a periodic system is the Floquet theory and it can be found in textbooks on ordinary differential equations [25,26]. An outline of the theory is given in the electronic supplementary material.

Let **x**(*t*) represents the vector of perturbations in the state variables[*δ*u *δ*w *δ*q *δ**θ*]^{T}. Equation (2.14) can be re-written as
2.15awhere **A**(*t*) is a periodic matrix with period 1, i.e.
2.15b

On the basis of the Floquet theory (see electronic supplementary material), the general solution of equation (2.15) can be written in the following form:
2.16where *a _{j}* is a constant,

**p**

*(*

_{j}*t*) (

*j*= 1, 2, 3, 4) are periodic vector functions with period 1 and

*ρ*

_{j}(

*j*= 1, 2, 3, 4) are the characteristic exponents of the system.

From equation (2.16), it is clear that when the real part of every *ρ*_{j} (*j* = 1, 2, 3, 4) is negative, **x**(*t*), the perturbation, will die out as time increases and the flight is stable; otherwise the flight is unstable; that is, the flight stability is dependent on *ρ*_{j} (*j* = 1, 2, 3, 4).

The characteristic exponents can be determined from **A**(*t*) in equation (2.15) as follows (see electronic supplementary material). First, solve the matrix equation
2.17aand
2.17bwhere ** ψ**(

*t*) is a matrix and each of its column vectors satisfies equation (2.15) (

**(**

*ψ**t*) is called as a fundamental matrix of equation (2.15) (see electronic supplementary material)) and

**I**is the unit matrix; equation (2.17) is solved numerically and matrix

**(**

*ψ**t*) is determined. Next, obtain the eigenvalues of

**(1),**

*ψ**ν*

_{j}(

*j*= 1, 2, 3, 4).

*ν*

_{j}can be written as

*ν*

_{j}= |

*ν*

_{j}|e

^{i}

*, where |*

^{φ}*ν*

_{j}| and

*φ*are the modulus and argument of ν

*, respectively. Finally, ρ*

_{j}*is determined as (see equation S21 of electronic supplementary material): 2.18*

_{j}## 3. Results and discussion

### 3.1. Periodic solution of hover flight

To obtain the periodic solution that represents the hover flight, the equations of motion coupled with the Navier–Stokes equations are solved numerically using the method by Wu *et al*. [8] outlined in electronic supplementary material. Quantitative assessment of the accuracy of the solution of the Navier–Stokes equations has been made to determine the grid size (see electronic supplementary material).

The numerical solution, *u*_{0}(*t*), *w*_{0}(*t*), *q*_{0}(*t*) and *θ*_{0}(*t*), for the hovering model dronefly is shown in figure 2*a*. As aforementioned, *α*_{d}, *α*_{u} and are determined in the solution process using the conditions of equilibrium flight: mean horizontal and vertical velocities of the body centre of mass and mean pitch rate of body are zero. It was found that for the hovering model dronefly, when *α*_{d} = 31.8°, *α*_{u} = 32.4° and , periodic solutions with zero means in *u*, *w* and *q* were obtained. The experimental measured values for *α*_{d}, *α*_{u} and are 33.8°, 33° and 7.12°, respectively [10]. The computed values are not very different from the experimental data. The solution for the hawkmoth is shown in figure 2*b*; the computed values of *α*_{d}, *α*_{u} and are 45.8°, 46.5° and 10.5°, respectively. The computed is not very different from the experimental value, 9.2° [21]. Using the computed *α*_{d} and *α*_{u}, the time-varying angle of attack is plotted in figure 3 (the *β* = 0 curve) compared with the experimental data; there is some discrepancy between the computed and experimental data (see below for more discussion on the discrepancy). It should be noted that in figure 2*b*, the mean of *w* is not exactly 0 but is about 0.005. That is, what we obtained is only an approximate solution of hovering. But the velocity is so small that the approximation is a good one.

For the hawkmoth, the peak-to-peak value of the horizontal velocity of the centre of mass of body is about 0.12*U* (figure 2*b*) and the peak-to-peak value of the pitch angle is about 5°. For the dronefly (figure 2*a*), the corresponding values are 0.025*U* and 0.3°. The body oscillations in the hawkmoth are much larger than that of the dronefly. This is expected because the wingbeat frequency of the hawkmoth is lower and wing mass (relative to the body mass) is much larger than that of the dronefly.

The periodic aerodynamic forces and moment [*XA*_{0}(*t*), *ZA*_{0}(*t*) and *MA*_{0}(*t*)] and wing-inertial forces and moment [*XI*_{0}(*t*) and *ZI*_{0}(*t*) and *MI*_{0}(*t*)] acting on the hovering insects are shown in figure 4*a* and *b* for the model dronefly and model hawkmoth, respectively. It is interesting to note that during the flight, the wing-inertial force and moment are of the same order of magnitude as the aerodynamic force and moment.

### 3.2. Aerodynamic and inertial derivatives and fundamental matrix

After the periodic solution of the hover flight is determined, aerodynamic and wing-inertial force and moment derivatives needed in matrix **A**(*t*) are computed. The aerodynamic derivatives with respect to *u* are computed as follows. Let the insect move at *u* = *u*_{0}(*t*) + 0.01, *w* = *w*_{0}(*t*), *q* = *q*_{0}(*t*) and *θ* = *θ*_{0}(*t*), and compute *XA*(*t*), *ZA*(*t*) and *MA*(*t*) and *XI*(*t*), *ZI*(*t*) and *MI*(*t*) at this motion; and then change *u* to *u* = *u*_{0}(*t*)−0.01 and compute the forces and moments. Based on these computed forces and moments, *δ**XA*(*t*)*δ**u*, *δ**ZA*(*t*)*δ**u*, etc., can be estimated using the finite-difference formulas (equation (2.13)). Derivatives with respect to *w* or *q* are computed similarly. The estimated derivatives are shown in electronic supplementary material figures S5 and S6. In the above computations, Δ*u* (or Δ*w* or Δ*q*) in equation (2.13) is taken as 0.02; we have tried other small values (0.04; 0.06) and the results were approximately the same.

With the aerodynamic and inertial derivatives computed, **A**(*t*) matrix of equation (2.14) or equation (2.15) is determined, and equation (2.17) is solved numerically to give the fundamental matrix, ** ψ**(

*t*). The matrix at

*t*= 1, which is used to compute the characteristic exponents of the system, is as follows.

For the model dronefly 3.1

For the model hawkmoth 3.2

### 3.3. Stability of the periodic solution

The eigenvalues of matrix ** ψ**(1) in equations (3.1) and (3.2) can be easily calculated. Putting the eigenvalues into equation (2.18) gives the characteristic exponents,

*ρ*

_{j}(

*j*= 1, 2, 3, 4), and they are shown in table 4. For each of the two model insects, there is a pair of complex characteristic exponents (

*ρ*

_{1,2}) with positive real part and two negative real characteristic exponents (

*ρ*

_{3}and

*ρ*

_{4}). Letting

*ρ*

_{1,2}=

*γ*± i

*ω*, the general solution of the equation of disturbance motion, equation (2.16), can be re-written as: 3.3

For both model insects, because *ρ*_{3} and *ρ*_{4} are negative, the second and third terms on the right-hand side (r.h.s.) of equation (3.3) die out as time increases; but because γ is positive, the first term on the r.h.s. of equation (3.3) will grow with time. Thus, the periodic solution, or the hover flight, is unstable.

### 3.4 Comparison with averaged-model theory

It is of interest to compare the above results of Floquet theory (cyclic-motion stability analysis) with those of the averaged-model theory (fixed-point stability analysis). A brief description of the averaged-model theory has been given in §1 and it is further described in the electronic supplementary material. In the averaged-model theory, the state variables are time-averaged over wingbeat cycle, denoted as , , and , and the equation of disturbance motion is a system of linear ordinary differential equations with constant coefficients. Techniques of eigenvalue and eigenvector analysis can be used to solve the equation.

The eigenvalues for the two hovering model insects have been calculated in the electronic supplementary material and the results are shown in table 5. For each of the insects, there are a pair of complex eigenvalues (*λ*_{1,2}) with a positive real part and two negative real eigenvalues (*λ*_{3} and *λ*_{4}). Let , the general solution of the equation of disturbance motion (equation S24) can be written as:
3.4
where *c _{j}* (

*j*= 1, 2, 3, 4) are constants and

**q**

*(*

_{j}*j*= 1, 2, 3, 4) are eigenvectors corresponding to

*λ*

_{j}(

*j*= 1, 2, 3, 4). The first term on the r.h.s. of equation (3.4) represents an unstable oscillatory divergence mode, the second term a stable fast convergence mode and the third term a stable slow convergence mode.

Comparing the results of the Floquet theory in the previous section with those of the averaged-model theory, the following observations can be made. For each of the two model insects, both theories predict the same unstable and stable natural modes of motion. For the model dronefly, the growth rate of the unstable mode and the decay rates of the stable modes given by the Floquet theory are very close to their counterparts of the averaged-model theory (*γ* = 0.0476, *ρ*_{3} = −0.115 and *ρ*_{4} = −0.0144 in the Floquet theory; , *λ*_{3} = −0.114 and *λ*_{4} = −0.0147 in the averaged-model theory). But for the model hawkmoth, the growth and decay rates given by the Floquet theory are significantly different from their counterparts of the averaged-model theory: for Floquet theory, *γ* = 0.199, *ρ*_{3} = −0.747 and *ρ*_{4} = −0.103, while for averaged-model theory, , *λ*_{3} = −0.715 and *λ*_{4} = −0.094. The real part of the root of the unstable mode given by the averaged-model theory ) is about 25 per cent larger than that by the Floquet theory (*γ* = 0.199).

### 3.5. Comparison with simulation using full equations of motion coupled with Navier–Stokes equations

In the present Floquet stability analysis, the Taylor expansions of the aerodynamic and wing-inertial forces and moments—equations (2.11) and (2.12)—contain only partial derivatives with respect to the state variables (*u*, *w* and *q*) and those with respect to , , , etc., are assumed to be negligible (this assumption is also made in the averaged-model theory). Here, we test the validity of the above assumption by comparing the results of the Floquet theory with those of simulation using the complete equations of motion coupled with Navier–Stokes equations. Results of the averaged-model theory are also included for comparison.

In the averaged-model theory, at equilibrium flight, , , and are zero. When initial disturbances are prescribed at *t*=0, the solution of the subsequent disturbance motion is given by equation (3.4). The eigenvectors in equation (3.4) (**q*** _{j}*,

*j*= 1, 2, 3 and 4) have been given in electronic supplementary material, table S3. The constants,

*c*

_{1},

*c*

_{2},

*c*

_{3}and

*c*

_{4}in the equation are determined as follows: let

*t*= 0 in equation (3.4); with , , and known, equation (3.4) becomes a linear algebraic equation, the solution of which gives

*c*

_{1},

*c*

_{2},

*c*

_{3}and

*c*

_{4}. Thus, the disturbance motion for the given initial values is determined analytically.

In the Floquet theory, at equilibrium flight, *u*_{0}(*t*), *w*_{0}(*t*), *q*_{0}(*t*) and *θ*_{0}(*t*) are known periodic functions of time. When the flight is disturbed, the subsequent disturbance motion is given by equation (3.3). But in this case, the vectors **p*** _{j}* (

*j*= 1,2,3,4) in equation (3.3) are not known (it is only known that they are periodic functions of time with the same period as that of the system matrix) and the constants

*a*(

_{j}*j*= 1,2,3,4) in equation (3.3) cannot be determined using the given initial values. That is, the Floquet theory can only give an analytical general solution (equation (3.3)) to the system equation (equation (2.15

*a*)) but cannot give an analytical specific solution at given initial conditions. However, by numerically solving equation (2.15

*a*), the specific solution at given initial conditions can be obtained. The initial conditions for equation (2.15

*a*) are: , , and ; i.e. we start from the periodic solution with the initial disturbances added to it.

In the simulation using the complete equations of motion (equation (2.1)) coupled with the Navier–Stokes equations (equation (2.2)), the same initial conditions as those used for equation (2.15*a*) in the Floquet theory are used.

First, we consider the model dronefly. In equation (3.4), if we choose the initial values , , and such that *c*_{3} = *c*_{4} = 0 and *c*_{1} and *c*_{2} are non-zero values, the disturbance motion would have only the first term, i.e. only the unstable divergence mode is excited; if we choose the initial values , , and such that *c*_{1} = *c*_{2} = *c*_{4} = 0 and *c*_{3} is a non-zero value, the disturbance motion would have only the term , i.e. only the fast stable subsidence mode is excited. Figure 5*a* shows the comparison of the results of the Floquet theory, the averaged-model theory and the numerical simulation when only the unstable divergence mode is excited. Figure 5*b* shows the corresponding results when only the fast stable subsidence mode is excited. Figure 5*c* shows the comparison of the results when the initial disturbances are , , which simulate a horizontal gust, and figure 5*d* shows the results when the initial disturbances are , , which simulate an initial pitch disturbance. From figure 5, we see that for the model dronefly, both the Floquet theory and the averaged-model theory are in good agreement with the numerical simulation.

Next, we consider the model hawkmoth. Computations with similar sets of initial conditions are conducted for the hawkmoth, and the results are shown in figure 6. It is seen that for the hawkmoth, the Floquet theory is in good agreement with the numerical simulation, while there are pronounced differences between the averaged-model theory and the simulation (figure 6*a*(iv), *c*, *d*(iv)).

Good agreement between the Floquet theory and the simulation for both the model insects indicates that the assumption made in the Floquet analysis (partial derivatives with respect to , , are negligible) is reasonable.

### 3.6. Model hawkmoth with tilted stroke plane

In the above, for comparison with dronefly (whose *β* is zero), *β* = 0 (and body angle *χ* = 54.8°) was used for the model hawkmoth. In fact, the hawkmoth considered hovers with *β* = 15° and *χ* = 39.8° [21]. We repeated the analysis of the model hawkmoth for the case of *β* = 15° (and *χ* = 39.8°).

The characteristic exponents of the Floquet theory are *ρ*_{1,2} = 0.181 ± 0.453i, *ρ*_{3} = −0.651 and *ρ*_{4} = −0.062 (the corresponding values for the case of *β* = 0 are 0.199 ± 0.562i, −0.747 and −0.103).

The eigenvalues of the averaged-model theory are *λ*_{1,2} = 0.210 ± 0.511i, *λ*_{3} = −0.641 and *λ*_{4} = −0.093 (the corresponding values for the case of *β* = 0 are 0.253 ± 0.588i, −0.715 and −0.094).

It is seen that when the stroke plane is tilted by 15°, the natural modes of motion are the same as those in the case of *β* = 0, but the growth of decay rate of a mode is decreased a little (e.g. in the Floquet analysis, the growth rate of the unstable mode is 0.181 for *β* = 15° and 0.199 for *β* = 0). But similar to the case of *β* = 0, for *β* = 15°, the growth and decay rates given by the Floquet theory are significantly different from their counterparts of the averaged-model theory; for example, the growth rate of the unstable mode given by the averaged-model theory (0.210) is about 16 per cent larger than that by the Floquet theory (0.181).

Figure 7 compares the results of the Floquet theory, the averaged-model theory and the simulation using full equations of motion coupled with the Navier–Stokes equations for the case of *β* = 15°. It is seen that similar to the case of *β* = 0, the Floquet theory is in good agreement with the numerical simulation, whilst there are relatively large differences between the averaged-model theory and the simulation.

### 3.7. Implications of the present results

#### 3.7.1. Cyclic-motion stability analysis needed for relatively large insects

From the above test of the theories by simulation using complete equations of motion coupled with the Navier–Stokes equations, we see that the Floquet theory (cyclic-motion stability analysis) agrees well with the simulation for both the model dronefly and the model hawkmoth, but the averaged-model theory (fixed-point stability analysis) gives good results only for the model dronefly. For the model hawkmoth, the Floquet theory (and the simulation) shows that in the unstable mode, *t*_{double} (time for an initial disturbance to double its value) is approximately 3.5 wingbeat cycles (0.693/0.199 ≈ 3.5); however, the value given by the averaged-model theory is 2.7 (0.693/0.253 ≈ 2.7), the growth rate of the instability being over-predicted by approximately 25 per cent by the fixed-point stability analysis. This shows that for a relatively large insect such as a hawkmoth, which has relatively low wingbeat frequency and large wing-mass to body-mass ratio, and hence pronounced cyclic-motion at wingbeat frequency, cyclic-motion stability analysis is required to give accurate prediction of its flight dynamic properties.

#### 3.7.2. Cyclic-motion equilibrium needs different control from that of fixed-point equilibrium

For a small insect like the model dronefly considered above, the averaged-model theory (fixed-point stability analysis) is applicable. In this case, when certain control force (or moment) is required, the force can be produced in any part of the wingbeat cycle, in the downstroke or in the upstroke, or a part of the force produced in the downstroke and the other part in the upstroke, as long as the mean force is equal to the required one.

But for an insect like the model hawkmoth, equilibrium flight involves pronounced cyclic motion and cyclic-motion stability analysis is required. In this case, the cyclic-motion at wingbeat frequency influences the free motion and the time-history of the control force in a wingbeat becomes important. For the control analysis of such insects, existing well-developed control theories for systems of fixed-point equilibrium are no longer applicable and new methods that take the cyclic variation of the flight dynamics into account are needed.

## Acknowledgements

This research was supported by grants from the National Natural Science Foundation of China (10732030) and the Foundation for Author of National Excellent Doctoral Dissertation of China (2007B31).

- Received January 29, 2012.
- Accepted March 13, 2012.

- This journal is © 2012 The Royal Society