The longitudinal dynamic flight stability of a hovering bumblebee was studied using the method of computational fluid dynamics to compute the aerodynamic derivatives and the techniques of eigenvalue and eigenvector analysis for solving the equations of motion.

For the longitudinal disturbed motion, three natural modes were identified:one unstable oscillatory mode, one stable fast subsidence mode and one stable slow subsidence mode. The unstable oscillatory mode consists of pitching and horizontal moving oscillations with negligible vertical motion. The period of the oscillations is 0.32 s (approx. 50 times the wingbeat period of the bumblebee). The oscillations double in amplitude in 0.1 s; coupling of nose-up pitching with forward horizontal motion (and nose-down pitching with backward horizontal motion) in this mode causes the instability. The stable fast subsidence mode consists of monotonic pitching and horizontal motions, which decay to half of the starting values in 0.024 s. The stable slow subsidence mode is mainly a monotonic descending (or ascending) motion, which decays to half of its starting value in 0.37 s.

Due to the unstable oscillatory mode, the hovering flight of the bumblebee is dynamically unstable. However, the instability might not be a great problem to a bumblebee that tries to stay hovering: the time for the initial disturbances to double (0.1 s) is more than 15 times the wingbeat period (6.4 ms), and the bumblebee has plenty of time to adjust its wing motion before the disturbances grow large.

In last twenty years, much work has been done to study the aerodynamics and energetics of insect flight, and considerable progress has been made in these areas (e.g. Dudley and Ellington, 1990a,b; Dickinson and Götz, 1993; Ellington et al., 1996; Dickinson et al., 1999; Wang, 2000; Sun and Tang, 2002; Usherwood and Ellington, 2002a,b). The area of insect flight stability has received much less consideration,however. Recently, with the current understanding of the aerodynamic force mechanisms of insect flapping wings, researchers are beginning to devote more effort to understanding this area.

Thomas and Taylor (2001)and Taylor and Thomas (2002)studied static stability (an initial directional tendency to return to equilibrium after a disturbance) of gliding animals and flapping flight,respectively. They found that flapping did not have any inherently destabilizing effect: beating the wing faster simply amplified the existing stability or instability, and that flapping could even enhance stability compared to gliding flight at a given speed.

Taylor and Thomas (2003)studied dynamic flight stability in the desert locust Schistocerca gregaria, providing the first formal quantitative analysis of dynamic stability in a flying animal (the dynamic stability of a flying body deals with the oscillation of the body about its equilibrium position following a disturbance). A very important assumption in their analysis was that the wingbeat frequency was much higher than the natural oscillatory modes of the insect, thus when analyzing its flight dynamics, the insect could be treated as a rigid flying body with only 6 degrees of freedom (termed rigid body approximation). In the rigid body approximation, the time variations of the wing forces and moments over the wingbeat cycle were assumed to average out;the effects of the flapping wings on the flight system were represented by the wingbeat-cycle average aerodynamic forces and moments that could vary with time over the time scale of the flying rigid body. In addition, the gyroscopic forces of the wings were assumed to be negligible. It was further assumed that the animal's motion consists of small disturbances from the equilibrium condition; as a result, the linear theory of aircraft flight dynamics was applicable to the analysis of insect flight dynamics. The authors first measured the aerodynamic force and moment variations of the tethered locust by varying the wind-tunnel speed and the attitude of the insect, obtaining the aerodynamic derivatives. Then they studied the longitudinal dynamic flight stability of the insect using the techniques of eigenvalue and eigenvector analysis.

In the study of Taylor and Thomas(2003), the dynamic stability of forward flight at high flight speed (the advance ratio was around 0.9) was studied. Many insects often hover. In hovering, unlike in forward flight, the stroke plane is generally approximately horizontal, the body angel is relatively large and the wing in the downstroke and in the upstroke operates under approximately the same conditions. As a result, the aerodynamic derivatives, hence the dynamic stability properties of hovering, must be different from those of forward flight. It is of great interest to investigate the dynamic flight stability of hovering.

In the present paper, we study the longitudinal dynamic flight stability in a hovering bumblebee. The bumblebee was chosen because previous studies on bumblebees provide the most complete morphological data and wing-motion descriptions. In the study by Taylor and Thomas(2003), due to the limits of the experimental conditions, the insect had to be tethered and the reference flight might not have been in equilibrium, so some derivatives could not be measured directly. If a computational method were used to obtain the aerodynamic derivatives, the above difficulties could be solved. More importantly, the computational approach allows simulation of the inherent stability of a flapping motion in the absence of active control. This is very difficult or impossible to achieve in experiments using real insects, as was done by Taylor and Thomas(2003). In the present study,we used the method of computational fluid dynamics (CFD) to compute the flows and to obtain the aerodynamic derivatives. First, conditions for force and moment equilibrium were determined. Then, the aerodynamic derivatives at equilibrium flight were computed. Finally, the longitudinal dynamic flight stability of the hovering bumblebee was studied using the techniques of eigenvalue and eigenvector analysis.

### Equations of motion

Similar to Taylor and Thomas(2003), we make the rigid body approximation: the insect is treated as a rigid body of 6 degrees of freedom(in the present case of symmetric longitudinal motion, only three degrees of freedom) and the action of the flapping wing is represented by the wingbeat-cycle average forces and moment (in addition, the gyroscopic effects of the wing are assumed negligible). This model of the hovering bumblebee is sketched in Fig. 1A.

Fig. 1.

(A) A sketch of the rigid body approximation. (B) Definition of the state variables u, w, q and θ. The bumblebee is shown during a perturbation (u, w, q and θ are zero at equilibrium).

Fig. 1.

(A) A sketch of the rigid body approximation. (B) Definition of the state variables u, w, q and θ. The bumblebee is shown during a perturbation (u, w, q and θ are zero at equilibrium).

Let oxyz be a non-inertial coordinate system fixed to the body. The origin o is at the center of mass of the insect and axes are aligned so that the x-axis is horizontal and points forward at equilibrium. The variables that define the motion(Fig. 1B) are the forward(u) and dorso-ventral (w) components of velocity along x- and z-axes, respectively, the pitching angular-velocity around the center of mass (q), and the pitch angle between the x-axis and the horizontal (θ). oExEyEzEis a coordinate system fixed on the earth; the xE-axis is horizontal and points forward.

The equations of motion are intrinsically non-linear, but may be linearized by approximating the body's motion as a series of small disturbance from a steady, symmetric reference flight condition. The linearized equations (see Etkin, 1972; Taylor and Thomas, 2003) are:
$\ {\delta}{\dot{u}}=X_{\mathrm{u}}{\delta}u{/}m+X_{\mathrm{w}}{\delta}w{/}m+X_{\mathrm{q}}{\delta}q{/}m-g{\delta}{\theta},$
1
$\ {\delta}{\dot{w}}=Z_{\mathrm{u}}{\delta}u{/}m+Z_{\mathrm{w}}{\delta}w{/}m+Z_{\mathrm{q}}{\delta}q{/}m,$
2
$\ {\delta}{\dot{q}}=M_{\mathrm{u}}{\delta}u{/}I_{\mathrm{y}}+M_{\mathrm{w}}{\delta}w{/}I_{\mathrm{y}}+M_{\mathrm{q}}{\delta}q{/}I_{\mathrm{y}},$
3
$\ {\delta}{\dot{{\theta}}}={\delta}q,$
4
$\ {\delta}{\dot{x}}_{\mathrm{E}}={\delta}u,$
5
$\ {\delta}{\dot{z}}_{\mathrm{E}}={\delta}w,$
6
where Xu, Xw, Xq, Zu, Zw, Zq, Mu, Mw and Mq are the aerodynamic derivatives [X and Z are the x- and z-components of the total aerodynamic force (due to the wing and the body), respectively, and M is the aerodynamic pitching moment (due to the wing and the body)]; m is the mass of the insect; g is the gravitational acceleration; Iy is the pitching moment of inertia about y axis; .' represents differentiation with respect to time (t); E and E represent the xE- and zE-component of the velocity of the mass center of the insect, respectively; the symbol δ denotes a small disturbance quantity. At reference flight (hovering), u, w, q, θ are zero (θ is zero because the x-axis is aligned with horizontal at reference flight), and X=0, Z=-mg and M=0(the forces and moments are in equilibrium).

In deriving the linearized equations (Equations 1, 2, 3, 4, 5, 6), the aerodynamic forces and moment (X, Z and M) are represented as analytical functions of the disturbed motion variables (δu, δw andδ q) and their derivatives(Etkin, 1972; Taylor and Thomas, 2003), e.g. X is represented as X=Xe+Xuδu+Xwδw+Xqδq,where the subscript e (for equilibrium) denotes the reference flight condition. In so doing, the effects of the whole body motion on the aerodynamic forces and moment are assumed to be quasisteady (terms that include δ, etc. are not included). The whole body motion is assumed to be slow enough for its unsteady effects to be negligible.

Let c, U and tw be the reference length,velocity and time, respectively [c is the mean chord length of the wing; U is the mean flapping velocity at the radius(r2) of the second moment of wing area, defined as U=2Φnr2 (Φ and n are the stroke amplitude and stroke frequency, respectively); tw is the period of the wingbeat cycle (tw= 1/n)]. The non-dimensional forms of Equations 1, 2, 3, 4, 5, 6 are:
$\ \left[\begin{array}{c}{\delta}{\dot{u}}^{+}\\{\delta}{\dot{w}}^{+}\\{\delta}{\dot{q}}^{+}\\{\delta}{\dot{{\theta}}}\end{array}\right]=\mathbf{A}\left[\begin{array}{c}{\delta}u^{+}\\{\delta}w^{+}\\{\delta}q^{+}\\{\delta}{\theta}\end{array}\right],$
7
$\ {\delta}{\dot{x}}_{\mathrm{E}}^{+}={\delta}u^{+},$
8
$\ {\delta}{\dot{z}}_{\mathrm{E}}^{+}={\delta}w^{+},$
9
where A is the system matrix:
$\ \mathbf{A}=\left[\begin{array}{cccc}X_{\mathrm{u}}^{+}a_{\mathrm{m}}{/}m&X_{\mathrm{w}}^{+}a_{\mathrm{m}}{/}m&X_{\mathrm{q}}^{+}a_{\mathrm{m}}{/}m&-g{/}a_{\mathrm{g}}\\Z_{\mathrm{u}}^{+}a_{\mathrm{m}}{/}m&Z_{\mathrm{w}}^{+}a_{\mathrm{m}}{/}m&Z_{\mathrm{q}}^{+}a_{\mathrm{m}}{/}m&0\\M_{\mathrm{u}}^{+}a_{\mathrm{I}}{/}I_{\mathrm{y}}&M_{\mathrm{w}}^{+}a_{\mathrm{I}}{/}I_{\mathrm{y}}&M_{\mathrm{q}}^{+}a_{\mathrm{I}}{/}I_{\mathrm{y}}&0\\0&0&1&0\end{array}\right],$
10
where am=0.5ρUSttw,
$$a_{\mathrm{I}}=0.5{\rho}U^{2}S_{\mathrm{t}}ct_{\mathrm{w}}^{2}$$
and ag = U/tw, and the non-dimensional forms are: δu+ = δu/Uw+ = δw/U, δq+ =δqtw, X =X/0.5ρU2St (ρ denotes the air density and St denotes the area of two wings),Z+=Z/0.5ρU2St, M+=Z/0.5ρU2Stc, t+=t/tw,
$${\delta}{\dot{x}}_{\mathrm{E}}^{+}={\delta}{\dot{x}}_{\mathrm{E}}{/}U$$
and
$${\delta}{\dot{z}}_{\mathrm{E}}^{+}={\delta}{\dot{z}}_{\mathrm{E}}{/}U$$
(using the flight data given below, am, aI and ag are computed as am=1.96 mg, aI=0.233×10-9 kg m2, ag=710.9 m s-1; ρ is 1.25 kg m-3and g is 9.8 m s-2).

### Flight data and non-dimensional parameters of wing motion

Flight data for the bumblebee are taken from Dudley and Ellington(1990a,b). The general morphological data are as follows: m=175 mg; wing length R=13.2 mm; c=4.01 mm, r2=0.55R;area of one wing (S) is 53 mm2; free body angle(χ0) is 57.5°; body length (lb) is 1.41R; distance from anterior tip of body to center of mass divided by body length (l) is 0.48lb, distance from wing base axis to center of mass divided by body length (ll) is 0.21 lb; pitching moment of inertia of the body about wing-root axis (Ib) is 0.48×10-8 kg m2. Assuming that the contribution of the wing mass to the pitching moment of inertia is negligible (the added-mass on the wings has been included in the CFD model), Iy, the pitching moment of the bumblebee about y-axis, can be computed as

$$I_{\mathrm{y}}=I_{\mathrm{b}}-l_{1}^{2}m^{2}=0.213{\times}10^{-8}{\\ }\mathrm{kg\ m}^{2}$$
⁠. Taylor and Thomas(2003) estimated the wings'contribution to the pitching moment of inertia for locusts and showed that the wings' contribution, which is proportional to the ratio of wing mass to the total body mass, was small, less than 3.5% of the pitching moment of inertia. The wings of locusts comprise around 4% of the total body mass; for the bumblebee, the wings comprise only 0.52% of the total body mass(Dudley and Ellington, 1990a). Therefore, the above estimation of Iy should be sufficient.

The wing-kinematic data are as follows: Φ=116°; n=155 Hz;β=6°; χ(body angle)=46.8°. U is computed as U=4.59 m s-1.

### Determination of the equilibrium conditions and computation of the aerodynamic derivatives

#### The wing, the body and the flapping motion

In determining the equilibrium conditions of the flight, we need to calculate the flows around the wings (at equilibrium, the body does not move);to obtain the aerodynamic derivatives, we need to compute the flows around the wing and around the body. In the present CFD model, it is assumed that the wings and body do not interact aerodynamically, neither do the contralateral wings, and the flows around the wings and body are computed separately. It is also assumed that the wing is inflexible. The wing planform used(Fig. 2A) is approximately the same as that of a bumblebee (Ellington,1984a). The wing section is a flat plate with rounded leading and trailing edges, the thickness of which is 3% of the mean chord length of the wing. The body of the insect is idealized as a body of revolution; the outline of the idealized body (Fig. 2B)is approximately the same as that of a bumblebee. Neglecting the axial asymmetry of the bumblebee can cause some differences in the computed body aerodynamic force. However, near hovering, the body aerodynamic force is much smaller than that of the wings (i.e. the aerodynamic force of the insect is mainly from the wings), and a small difference in the body aerodynamic force may not affect the aerodynamic derivatives greatly (see below).

Fig. 2.

(A) Portions of the grid for the bumblebee wing. (B) The bumblebee body planform and a portion of the grid.

Fig. 2.

(A) Portions of the grid for the bumblebee wing. (B) The bumblebee body planform and a portion of the grid.

The flapping motion of the wing (Fig. 3) consists of two parts: the translation (azimuthal rotation) and the rotation (flip rotation, or rotation around an axis along the wing). The velocity at the span location r2 due to wing translation is called the translational velocity (ut). The azimuth-rotational velocity of the wing(
$${\dot{{\phi}}}$$
) is related to
$$u_{\mathrm{t}}:{\\ }{\dot{{\phi}}}({\tau})=u_{\mathrm{t}}{/}r_{2}$$
,where τ is non-dimensional time. For the bumblebee, on the basis of data given by Dudley and Ellington(1990a), utis approximated by the simple harmonic function:
$\ u_{\mathrm{t}}^{+}=0.5{\pi}\mathrm{sin}(2{\pi}{\pi}{/}{\pi}_{\mathrm{c}}),$
11
where the non-dimensional translational velocity
$$u_{\mathrm{t}}^{+}=u_{\mathrm{t}}{/}U$$
,non-dimensional time τ=tU/c, and τc is the non-dimensional period of a wingbeat cycle. The geometric angle of attack of a wing is denoted by α. On the basis of flight data(Dudley and Ellington, 1990a; Ellington, 1984b), time variations of α are approximated as follows. α takes a constant value except at the beginning or near the end of a half-stroke. The constant value, called as mid-stroke angle of attack, is denoted by αdfor the downstroke and αu for the upstroke. Around stroke reversal, the wing flips and α changes with time. The angular velocity(
$${\dot{{\alpha}}}$$
) is given by:
$\ \begin{array}{c}{\dot{{\alpha}}}^{+}=0.5{\dot{{\alpha}}}_{0}^{+}\{1-\mathrm{cos}[2{\pi}({\tau}-{\tau}_{\mathrm{r}}){/}{\Delta}{\tau}_{\mathrm{r}}]\};\\{\tau}_{\mathrm{r}}{\leq}{\tau}{\leq}({\tau}_{\mathrm{r}}+{\Delta}{\tau}_{\mathrm{r}}),\end{array}$
12
where the non-dimensional form
$${\dot{{\alpha}}}^{+}={\dot{{\alpha}}}c{/}U$$
,
$${\dot{{\alpha}}}_{0}^{+}$$
is a constant, τr is the non-dimensional time at which the flip rotation starts, and Δτr the non-dimensional time interval over which the flip rotation lasts. In the time interval ofΔτ r, the wing rotates from α=αd toα=180°-αu. Therefore, when αdu and Δτr are specified,
$${\dot{{\alpha}}}_{0}^{+}$$
can be determined (around the next stroke reversal, the wing would rotate fromα=180°-αu to α=αd, the sign of the right-hand side of Equation 12 should be reversed). Δτr is termed flip duration. It is assumed that the axis of the pitching rotation is located at 0.3c from the leading edge of the wing.
Fig. 3.

Sketches of the reference frames and wing motion. o1x1y1z1is an inertial frame, with the x1,y1plane in the horizontal plane. oxyz′ is another inertial frame, with the xy′ plane in the stroke plane. φ, positional angle of the wing; φminand φmax, minimum and maximum positional angle, respectively;α, angle of attack of the wing; β, stroke plane angle; R,wing length. CL,w and CT,w,coefficients of vertical and thrust of wing, respectively.

Fig. 3.

Sketches of the reference frames and wing motion. o1x1y1z1is an inertial frame, with the x1,y1plane in the horizontal plane. oxyz′ is another inertial frame, with the xy′ plane in the stroke plane. φ, positional angle of the wing; φminand φmax, minimum and maximum positional angle, respectively;α, angle of attack of the wing; β, stroke plane angle; R,wing length. CL,w and CT,w,coefficients of vertical and thrust of wing, respectively.

In the flapping motion described above, the mid-stroke angles of attack(αd and αu), the flip duration(Δτr), the flip timing (τr), the period of flapping cycle (τc), the mean positional angle(

$${\bar{{\phi}}}$$
⁠) and the stroke plane angle(β) must be given. The Reynolds number (Re), which appears in the non-dimensional Navier-Stokes equations, is defined as Re=Uc/ν (ν is the kinematic viscosity of the air). The non-dimensional kinematic parameters are computed below.

U has been computed above. Re and τc are computed as Re=1326 and τc=7.12. On the basis of the flight data in Dudley and Ellington(1990a), the flip duration(Δτr) is set to 0.22τc and the flip rotation is assumed to be symmetrical (thus the flip timing τris determined in terms of Δτr). αdu and

$${\bar{{\phi}}}$$
are yet to be specified.

#### The flow solution method and evaluation of aerodynamic forces and moments

The flow equations (the Navier-Stokes equations) and the solution method used in the present study are the same as those described in Sun and Tang(2002). Once the flow equations are numerically solved, the fluid velocity components and pressure at discretized grid points for each time step are available. The aerodynamic forces and moments acting on the wing (or the body) are calculated from the pressure and the viscous stress on the wing (or the body) surface.

Resolving resultant aerodynamic force of the wing into the z1- and x1-axes, we obtain the vertical (Lw) and the horizontal (Tw)forces of the wing, respectively (see Fig. 3B,C; note that when computing aerodynamic derivatives with respect to q, the stroke plane and x1- and z1-axes rotate about the center of mass of the insect, Lw is not in vertical direction and Twnot in horizontal direction). Let my1,w be the moment about the y1-axis (which passes the wing root). The pitching moment about the center of mass of the insect due to the aerodynamic force of the wing (my,w) can be calculated (see Ellington, 1984a) as:
$\ M_{\mathrm{y,w}}=m_{\mathrm{y}_{1},\mathrm{w}}+L_{w}{\cdot}l_{1}{\\ }\mathrm{sin}({\chi}_{0}-{\chi})+T_{w}{\cdot}l_{1}{\\ }\mathrm{cos}({\chi}_{0}-{\chi}).$
13

The lift (Lb) and drag (Db) of the body are the vertical (z1 direction) and horizontal(x1 direction) components of the resultant aerodynamic force of the body, respectively. The pitching moment of the body(my,b) is the moment about the mass center due to the aerodynamic force of the body. The above forces and moment are non-dimensionalized by 0.5ρU2Stand 0.5ρU2Stc,respectively. The coefficients of Lw, Tw, my,w, Lb, Db and my,b are denoted a CL,w, CT,w, CM,w, CL,b, CD,b and CM,b, respectively.

#### Force and moment equilibrium

As seen above, the kinematic parameters of the wing left undetermined are the mid-stroke angles of attack (αd, αu) and the mean positional angle of the wing(

$${\bar{{\phi}}}$$
⁠). In the present study,α d, αu and
$${\bar{{\phi}}}$$
are not treated as known input parameters but are determined in the calculation by the force balance and moment balance conditions, i.e. the mean vertical force of the wings is equal to insect weight, the mean horizontal force of the wings is equal to zero, and the mean pitching moment of the wings (about the mass center) is equal to zero. The non-dimensional weight of the insect is defined as mg/0.5ρU2St, and its value is computed as 1.25. The mean vertical force coefficient of the wing needs to equal 1.25.

#### Aerodynamic derivatives

Conditions in equilibrium flight are taken as the reference conditions in the aerodynamic derivative calculations. In order to estimate the partial derivatives Xu, Xw, Xq, Zu, Zw, Zq, Mu, Mw and Mq, we make three consecutive flow computations for the wing: a u-series, in which u+ is varied whilst w+, q+ and θ are fixed at the reference values (i.e. w+, q+ andθ are zero), a w-series, in which w+ is varied whilst w+, q+ and θ are fixed at zero and a q-series, in which w+, q+ and θ are fixed at zero (in all the three series,wing kinematical parameters are fixed at the reference values); similar flow computations are conducted for the body. Using the computed data, curves representing the variation of the aerodynamic forces and moments with each of the w+, q+ and θ variables are fitted. The partial derivatives are then estimated by taking the local tangent(at equilibrium) of the fitted curves.

### Solution of the small disturbance equations

After the aerodynamic derivatives are determined, the elements of the system matrix A would be known. Equation 7 can be solved to yield insights into the dynamic flight stability of the hovering bumblebee.

The general theory for such a system of linear equations is in any textbook on flight dynamics (e.g. Etkin,1972); a concise description of the theory can be found in Taylor and Thomas (2003). Only an outline of the theory is given here. The central elements of the solutions for free motion, i.e. of the dynamic stability problem, are the eigenvalues and eigenvectors of A. In general, a l×l real matrix has l eigenvalues (λ12,... λ1) and l corresponding eigenvectors; an eigenvalue can be a real number (the corresponding eigenvector is real) or a complex number (the corresponding eigenvector is complex), and the complex eigenvalues (and eigenvectors) occur in conjugate pairs. A real eigenvalue and the corresponding eigenvector (or a conjugate pair of complex eigenvalues and the corresponding eigenvector pair) represent a simple motion called natural mode of the system. The free motion of the flying body after an initial deviation from its reference flight is a linear combination of the natural modes. Therefore, to know the dynamic stability properties of the system, one only needs to examine the motions represented by the natural modes. In a natural mode, the real part of the eigenvalue determines the time rate of growth of the disturbance quantities and the eigenvector determines the magnitudes and phases of the disturbance quantities relative to each other. A positive real eigenvalue will result in exponential growth of each of the disturbance quantities, so the corresponding natural mode is dynamically unstable (termed unstable divergent mode). The time to double the starting value is given by:
$\ t_{\mathrm{double}}=0.693{/}{\lambda}{\\ }({\lambda}{>}0).$
14
A negative eigenvalue will result in exponential decay of the disturbance quantities and the corresponding natural mode is dynamical stable (termed stable subsidence mode). The time to half the starting value is given by:
$\ t_{\mathrm{half}}=0.693\|{\lambda}|{\\ }({\lambda}{<}0).$
15
A pair of complex conjugate eigenvalues, e.g.
$${\lambda}_{1,2}={\hat{n}}{\pm}{\hat{{\omega}}}i$$
,will result in oscillatory time variation of the disturbance quantities with
$${\hat{{\omega}}}$$
as its angular frequency; the motion decays when is negative(dynamical stable; termed stable oscillatory mode) but grows when is positive (dynamical unstable;termed unstable oscillatory mode). The period (T) of the oscillatory motion is:
$\ T=2{\pi}{/}{\hat{{\omega}}},$
16
and the times to double or half the oscillatory amplitude are
$\ t_{\mathrm{double}}{\\ }\mathrm{or}{\\ }t_{\mathrm{half}}=0.693{/}|{\hat{n}}|.$
17
See fig. 10 of Taylor and Thomas(2003) for sketches of the four types of solution to the small disturbance equations.

The solution process of the present problem is summarized as follows. The eigenvalues and eigenvectors of A in Equation 7 are calculated, giving the natural modes; analyzing the motions of the natural modes gives the dynamic stability properties of the hovering bumblebee.

### Code validation and grid resolution test

The code used for the flow computations is the same as that in Sun and Tang(2002) and Sun and Wu(2003). It was tested in Sun and Wu (2003) using measured unsteady aerodynamic forces on a flapping model fruitfly wing. The calculated drag coefficient agreed well with the measured value. For the lift coefficient, the computed value agreed well with the measured value, except at the beginning of a half-stroke, where the computed peak value was smaller than the measured value. The discrepancy might be because the CFD code does not resolve the complex flow at stroke reversal satisfactorily. There is also the possibility that it is due to variations in the precise kinematic patterns,especially at the stroke reversal. Wu and Sun(2004) further tested the code using the recent experimental data by Usherwood and Ellington(2002b) on a revolving model bumblebee wing. In the whole α range (from -20° to 100°), the computed lift coefficient agreed well with the measured values. The computed drag coefficient also agreed well with the measured values except when αis larger than approximately 60°.

In the above computations, the computational grid was of the O-H type and had 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. The outer boundary was set at 20 chord lengths from the wing. The time step was 0.02. 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.

In the present study for the wing, we used similar grid dimensions as used in the test calculations (Wu and Sun,2004); for the body, the grid dimensions were 71×73×96 in the normal direction, along the body axis and in the azimuthal direction,respectively (tests have been conducted to show that these grid dimensions are appropriate for the present computations).

### The equilibrium flight

For different set of values of αd, αu and

$${\bar{{\phi}}}$$
⁠, the mean vertical and horizontal forces and mean pitching moment of the wings would be different.α d, αu and
$${\bar{{\phi}}}$$
are determined using equilibrium conditions. The calculation proceeds as follows. A set of values forα d, αu and
$${\bar{{\phi}}}$$
is guessed; the flow equations are solved and the corresponding mean vertical force(L,w), mean horizontal force (T,w) and mean moment (M,w) coefficients of the wing are calculated. If L,w is not equal to 1.25(the non-dimensional weight), or T,w is not equal to zero,or M,w is not equal to zero, αd, αu and
$${\bar{{\phi}}}$$
are adjusted; the calculations are repeated until the magnitudes of difference between L,w and 1.25, between T,w and 0 and between M,w and 0 are less than 0.01. The calculated results show that when, αd=27°,α u=21° and
$${\bar{{\phi}}}$$
=1°, the equilibrium conditions are satisfied.

### The aerodynamic derivatives

As defined above, X+ and Z+ are x- and z-components of the non-dimensional total aerodynamic force due to the wing and the body and M+ is the corresponding non-dimensional pitching moment. After the equilibrium flight conditions have been determined, aerodynamic forces and moments on the wing and on the body for each of u, w and q varying independently from the equilibrium value are computed. The corresponding X+, Z+ and M+ are obtained. In Fig. 4A-C, the u-series, w-series and q-series data, respectively,are plotted (in the figure the equilibrium value has been subtracted from each quantity). X+, Z+ and M+ vary approximately linearly with u+w+ and q+ in a range of-0.1≤Δu+, Δw+ andΔ q+≤0.1, showing that the linearization of the equations of motion is only justified for small disturbances. (In the computations, we found that the aerodynamic forces and moment of the body are negligibly small compared to those of the wing; this is because the relative velocity that the body sees is very small.) The aerodynamic derivatives

$$X_{\mathrm{u}}^{+}$$
⁠,
$$Z_{\mathrm{u}}^{+}$$
,
$$M_{\mathrm{u}}^{+}$$
,
$$X_{\mathrm{w}}^{+}$$
,
$$Z_{\mathrm{w}}^{+}$$
,
$$M_{\mathrm{w}}^{+}$$
,
$$X_{\mathrm{q}}^{+}$$
,
$$Z_{\mathrm{q}}^{+}$$
and
$$M_{\mathrm{q}}^{+}$$
, estimated using the data in Fig. 4, are shown in Table 1.

Fig. 4.

The u-series (A), w-series (B) and q-series (C)force and moment data. ΔX+ andΔ Z+, non-dimensional x- and z-components of the total aerodynamic force, respectively;Δ M+, non-dimensional pitching moment;Δ u+ and Δw+,non-dimensional x- and z-components of velocity of center of mass, respectively; Δq+, non-dimensional pitching rate (the equilibrium value is subtracted from each quantity).

Fig. 4.

The u-series (A), w-series (B) and q-series (C)force and moment data. ΔX+ andΔ Z+, non-dimensional x- and z-components of the total aerodynamic force, respectively;Δ M+, non-dimensional pitching moment;Δ u+ and Δw+,non-dimensional x- and z-components of velocity of center of mass, respectively; Δq+, non-dimensional pitching rate (the equilibrium value is subtracted from each quantity).

Table 1.

Non-dimensional aerodynamic derivatives

$$X_{\mathrm{u}}^{+}$$
$$Z_{\mathrm{u}}^{+}$$
$$M_{\mathrm{u}}^{+}$$
$$X_{\mathrm{w}}^{+}$$
$$Z_{\mathrm{w}}^{+}$$
$$M_{\mathrm{w}}^{+}$$
$$X_{\mathrm{q}}^{+}$$
$$Z_{\mathrm{q}}^{+}$$
$$M_{\mathrm{q}}^{+}$$
-0.785 -0.031 2.389 0.050 -1.033 -0.190 -0.090 -0.031 -0.883
$$X_{\mathrm{u}}^{+}$$
$$Z_{\mathrm{u}}^{+}$$
$$M_{\mathrm{u}}^{+}$$
$$X_{\mathrm{w}}^{+}$$
$$Z_{\mathrm{w}}^{+}$$
$$M_{\mathrm{w}}^{+}$$
$$X_{\mathrm{q}}^{+}$$
$$Z_{\mathrm{q}}^{+}$$
$$M_{\mathrm{q}}^{+}$$
-0.785 -0.031 2.389 0.050 -1.033 -0.190 -0.090 -0.031 -0.883

$$X_{\mathrm{u}}^{+},{\\ }Z_{\mathrm{u}}^{+}$$
and Mu+, non-dimensional derivatives of the x- and z-component of the aerodynamic force and aerodynamic moment, respectively, with respect to the x-component(u+) of the non-dimensional velocity;
$$X_{\mathrm{w}}^{+},{\\ }Z_{\mathrm{w}}^{+}$$
and
$$M_{\mathrm{w}}^{+}$$
, non-dimensional derivatives of the x- and z-component of the aerodynamic force and aerodynamic moment, respectively, with respect to the z-component(w+) of the non-dimensional velocity;
$$X_{\mathrm{q}}^{+},{\\ }Z_{\mathrm{q}}^{+}$$
and
$$M_{\mathrm{q}}^{+}$$
, non-dimensional derivatives of the x- and z-component of the aerodynamic force and aerodynamic moment, respectively, with respect to the non-dimensional pitching rate (q+).

Let us examine the aerodynamic derivatives and discuss how they are produced. First, we consider the derivatives with respect to u+. As seen in Table 1,

$$Z_{\mathrm{u}}^{+}$$
is almost zero,
$$X_{\mathrm{u}}^{+}$$
is negative and
$$M_{\mathrm{u}}^{+}$$
is positive and large. Fig. 5 shows the differences between CL,w, CT,w and CM,w at Δu+=0.05(Δu+q+=0) and their counterparts at reference flight. For convenience, we define a non-dimensional time, τ, such that τ=0 at the start of the downstroke and τ=1 at the end of subsequent upstroke. Differences between CL,w, CT,w and CM,w in some flight conditions and their counterparts at reference flight are denoted as ΔCL,wCT,w and ΔCM,w,respectively. In the reference flight (hovering), the stroke plane is almost horizontal. When the insect moves forward with Δu+,in the downstroke, the wing sees a larger relative velocity than that in the reference flight and its drag is larger than the reference value(ΔCT,w positive, see Fig. 5B), resulting in a decrease in X+; in the upstroke, the wing sees a smaller velocity than that in the reference flight and its drag is smaller than the reference value (ΔCT,w also positive, Fig. 5B), also resulting a decrease in X+. This explains the negative
$$X_{\mathrm{u}}^{+}$$
. As for the vertical force, there is an increase in the downstroke and a decrease in the upstroke compared to the reference value (see Fig. 5A),resulting in little change in Z+, which explains the small
$$Z_{\mathrm{u}}^{+}$$
. Since the wing is above the mass center, the decrease in X+ produces a nose-up pitching moment. From Fig. 5A, it is seen that ΔCL,w in the second half of the downstroke(τ=0.25-0.5) is larger than ΔCL,w in the first half of the downstroke (τ=0-0.25), producing a couple-nose-up pitching moment; similarly, ΔCL,w in the upstroke produces also produces a nose-up pitching moment. This explains the large positive
$$M_{\mathrm{u}}^{+}$$
.

Fig. 5.

Time courses of ΔCL,w (A),Δ CT,w (B) and ΔCM,w (C) in one flapping cycle. ΔCL,wCT,w and ΔCM,w are the difference between CL,w, CT,w and CM,w at Δu+=0.05(Δw+q+=0) and their counterparts at reference flight, respectively. Δu+and Δw+, non-dimensional x- and z-components of velocity of center of mass, respectively;Δ q+, non- dimensional pitching rate; τ,non-dimensional time.

Fig. 5.

Time courses of ΔCL,w (A),Δ CT,w (B) and ΔCM,w (C) in one flapping cycle. ΔCL,wCT,w and ΔCM,w are the difference between CL,w, CT,w and CM,w at Δu+=0.05(Δw+q+=0) and their counterparts at reference flight, respectively. Δu+and Δw+, non-dimensional x- and z-components of velocity of center of mass, respectively;Δ q+, non- dimensional pitching rate; τ,non-dimensional time.

Next, we examine the derivatives with respect to w+.

$$X_{\mathrm{w}}^{+}$$
and
$$M_{\mathrm{w}}^{+}$$
are very small and
$$Z_{\mathrm{w}}^{+}$$
is relatively large(Table 1). Fig. 6 shows the differences between CL,w, CT,w and CM,w at Δw+=0.05(Δu+q+=0) and their counterparts at reference flight. When the insect moves with positiveΔ w+, the wing sees an upward velocity in both the down- and upstrokes, and the lift, drag and moment of the wing in both the half-strokes would be increased compared to those at reference flight. As a result, ΔCL,w in both the down- and upstrokes is positive (Fig. 6A), resulting in a relatively large decrease in
$$Z_{\mathrm{w}}^{+}$$
CT,w in the downstroke has different sign from that in the upstroke (Fig. 6B) and so does ΔCM,w(Fig. 6C), resulting in small values in
$$X_{\mathrm{w}}^{+}$$
and
$$M_{\mathrm{w}}^{+}$$
. This explains the relatively large
$$Z_{\mathrm{w}}^{+}$$
and small
$$X_{\mathrm{w}}^{+}$$
and
$$M_{\mathrm{w}}^{+}$$
.

Fig. 6.

Time courses of ΔCL,w (A),Δ CT,w (B) and ΔCM,w (C) in one flapping cycle. ΔCL,wCT,w and ΔCM,w are the difference between CL,w, CT,w and CM,w at Δw+=0.05(Δu+q+=0) and their counterparts at reference flight, respectively. Δu+and Δw+, non-dimensional x- and z-components of velocity of center of mass, respectively;Δ q+, non-dimensional pitching rate; τ,non-dimensional time.

Fig. 6.

Time courses of ΔCL,w (A),Δ CT,w (B) and ΔCM,w (C) in one flapping cycle. ΔCL,wCT,w and ΔCM,w are the difference between CL,w, CT,w and CM,w at Δw+=0.05(Δu+q+=0) and their counterparts at reference flight, respectively. Δu+and Δw+, non-dimensional x- and z-components of velocity of center of mass, respectively;Δ q+, non-dimensional pitching rate; τ,non-dimensional time.

Finally, we examine the derivatives with respect to q+. As seen in Table 1,

$$X_{\mathrm{q}}^{+}$$
and
$$Z_{\mathrm{q}}^{+}$$
are very small and
$$M_{\mathrm{q}}^{+}$$
is relatively large. Fig. 7 shows the differences between CL,w, CT,w and CM,w at Δq+=0.07(Δu+w+=0) and their counterparts at reference flights. ΔCL,w andΔ CT,w are very small everywhere in the downstroke and the upstroke (Fig. 7A,B),resulting in the very small
$$X_{\mathrm{q}}^{+}$$
and
$$Z_{\mathrm{q}}^{+}$$
. ΔCM,w is small in a large part of the downstroke or the upstroke but is relatively large near the end of the half-stroke (Fig. 7C), resulting in the relatively large
$$M_{\mathrm{q}}^{+}$$
. Note that near the end of the half-strokes, ΔCL,w andΔ CT,w are very small, whereasΔ CM,w is relatively large. This means that the position of the action-line of the aerodynamic force of the wing must be changed by the whole body rotation of the insect (i.e. by the rotation of the stroke plane) near the end of the half-strokes, when the wing is in flip rotation. We thus see that although the whole body rotation of the insect could not change the magnitude of the aerodynamic force of the wing greatly from that of the reference flight, it changes the position of the action-line of the aerodynamics force near the end of the half-strokes, producing a pitching moment.

Fig. 7.

Time courses of ΔCL,w (A),Δ CT,w (B) and ΔCM,w (C) in one flapping cycle. ΔCL,wCT,w and ΔCM,w are the difference between CL,w, CT,w and CM,w at Δq+=0.07(Δu+w+=0)) and their counterparts at reference flight, respectively. Δu+and Δw+, non-dimensional x- and z-components of velocity of center of mass, respectively;Δ q+, non- dimensional pitching rate; τ,non-dimensional time.

Fig. 7.

Time courses of ΔCL,w (A),Δ CT,w (B) and ΔCM,w (C) in one flapping cycle. ΔCL,wCT,w and ΔCM,w are the difference between CL,w, CT,w and CM,w at Δq+=0.07(Δu+w+=0)) and their counterparts at reference flight, respectively. Δu+and Δw+, non-dimensional x- and z-components of velocity of center of mass, respectively;Δ q+, non- dimensional pitching rate; τ,non-dimensional time.

### The eigenvalues and eigenvectors

With the aerodynamic derivatives computed, the elements in the system matrix A are now known. The eigenvalues and the corresponding eigenvectors can then be computed, and the results are shown in Tables 2 and 3.

Table 2.

Eigenvalues of the system matrix

Mode 1, λ1,2Mode 2, λ3Mode 3, λ4
0.045±0.129i -0.197 -0.012
Mode 1, λ1,2Mode 2, λ3Mode 3, λ4
0.045±0.129i -0.197 -0.012

λ1,2, a pair of complex conjugate eigenvalues.λ 3 and λ4, real eigenvalues.

i, imaginary number.

Table 3.

Eigenvectors of the system matrix

Mode 1Mode 2Mode 3
δu+ 0.035±0.091i -0.071 -0.079
δw+ -4.95×10-4±2.04×10-4i 2.29×10-4 0.995
δq+ 0.121±0.061i 0.193 -6.57×10-4
δθ 0.706±0.689i -0.979 0.057
Mode 1Mode 2Mode 3
δu+ 0.035±0.091i -0.071 -0.079
δw+ -4.95×10-4±2.04×10-4i 2.29×10-4 0.995
δq+ 0.121±0.061i 0.193 -6.57×10-4
δθ 0.706±0.689i -0.979 0.057

δu+, δw+q+ and δθ, disturbance quantities in non-dimensional x-component and z-component of velocity,pitching rate and pitch angle, respectively. i, imaginary number.

As seen in Table 2, there are a pair of complex eigenvalues with a positive real part and two negative real eigenvalues, representing an unstable oscillatory motion (mode 1) and two stable subsidence motions (mode 2 and mode 3), respectively. The period for the oscillatory mode and the tdouble or thalf for the three modes, computed using Equations 14, 15, 16, 17, are shown in Table 4. Hereafter, we call modes 1, 2 and 3 unstable oscillatory mode, fast subsidence mode and slow subsidence mode, respectively.

Table 4.

Non-dimensional time constants of the natural modes

Mode 1
Mode 2
Mode 3
StabilityTtdoubleStabilitythalfStabilitythalf
Unstable 48.7 15.4 Stable 3.5 Stable 57.8
Mode 1
Mode 2
Mode 3
StabilityTtdoubleStabilitythalfStabilitythalf
Unstable 48.7 15.4 Stable 3.5 Stable 57.8

T, non-dimensional period of the oscillatory mode; tdouble, non-dimensional time for a divergence motion to double in amplitude; thalf, non-dimensional time for a damped motion to halve in amplitude.

As mentioned in Materials and methods, the eigenvector determines the magnitudes and phases of the disturbance quantities relative to each other. These properties can be clearly displayed by expressing the eigenvector in polar form (Table 5): since the actual magnitude of an eigenvector is arbitrary, only its direction is unique,and we have scaled them to make δθ=1.

Table 5.

Magnitudes and phase angles of the components of each of the three eigenvectors

Modeδu+δw+δq+δθ
Unstable oscillatory 1.0×10-1 (113.3°) 1.0×10-3 (156.7°) 1.4×10-1 (71.1°) 1 (0°)
Fast subsidence 0.7×10-1 (0°) 2.4×10-4 (180°) 2.0×10-1 (180°) 1 (0°)
Slow subsidence 1.4 (0°) 1.7×10 (0°) 1.2×10-2 (180°) 1 (0°)
Modeδu+δw+δq+δθ
Unstable oscillatory 1.0×10-1 (113.3°) 1.0×10-3 (156.7°) 1.4×10-1 (71.1°) 1 (0°)
Fast subsidence 0.7×10-1 (0°) 2.4×10-4 (180°) 2.0×10-1 (180°) 1 (0°)
Slow subsidence 1.4 (0°) 1.7×10 (0°) 1.2×10-2 (180°) 1 (0°)

δu+, δw+q+ and δθ, disturbance quantities in non-dimensional x-component and z-component of velocity,pitching rate and pitch angle, respectively.

Numbers in the parentheses are phase angles.

### The unstable oscillatory mode

The non-dimensional period of the oscillatory mode is T=48.7 and the non-dimensional time of doubling the amplitude is Tdouble=15.4 (Table 4). Note that the reference time used in non-dimensionalization of the equations of motion is the period of wingbeat cycle. Thus the period of the insect oscillation is about 49 times of the wingbeat period (the wingbeat period is tw=1/n=6.4 ms), and the starting value of the oscillation will double in 15 wingbeats.

As seen in Table 5, the unstable oscillatory mode is a motion in which δq+,δθ and δu+ are the main variables(δw+ is smaller than δq+and δu+ by two orders of magnitude; δθis seen to be very large, but it is the result ofδ q+ and the long period).δ q+ represents pitching motion and

$${\delta}u^{+}={\delta}{\dot{x}}_{\mathrm{E}}^{+}$$
(Equation 8) represents horizontal motion. Thus, in this mode the bumblebee conducts horizontal and pitching oscillations (it should be pointed out that in general, the motion of the system is a linear superposition of the simple motions represented by the natural modes, but if the initial conditions are correctly chosen, the motion represented by a natural mode can occur). The characteristic transients ofδ u+(
$${\delta}{\dot{x}}_{\mathrm{E}}^{+}$$
),δ q+ and δθ in this mode are plotted in Fig. 8A. It is seen that in a large part of a cycle, the bumblebee pitches down while moving backwards or pitches up while moving forward. The motion is sketched in Fig. 8B. As discussed below,pitching down (or up) while moving backwards (or forward) has a large destabilizing effect.

Fig. 8.

(A) Characteristic transients of disturbance quantities in the unstable oscillatory mode; (B) sketches showing the combinations of the pitching and horizontal motions in this mode. δu+q+ and δθ, disturbance quantities in non-dimensional x-component of velocity, pitching rate and pitch angle, respectively.

Fig. 8.

(A) Characteristic transients of disturbance quantities in the unstable oscillatory mode; (B) sketches showing the combinations of the pitching and horizontal motions in this mode. δu+q+ and δθ, disturbance quantities in non-dimensional x-component of velocity, pitching rate and pitch angle, respectively.

### The fast subsidence mode

For the fast subsidence mode, thalf is 3.5(Table 4); disturbances decrease to half of the starting values in about four wingbeats. As seen in Table 5, in this modeδ q+, θ and δu+ are also the main variables (δw+ is smaller by 3 orders of magnitude). δq+ and θ are out of phase(they have opposite signs); δu+ and θ are in phase. That is, when δθ has a positive initial value, so doesδ u+, but δq+ has a negative initial value. The insect would pitch down (back to the reference attitude) and at the same time moves forward (see the sketch in Fig. 9B). Note that this is different from the case of the unstable oscillatory mode, in which the insect pitches down while moving backwards and pitches up while moving forward (see Fig. 8B). The characteristic transients of δu+, δq+ andθ are plotted in Fig. 9A.

Fig. 9.

(A) Characteristic transients of disturbance quantities in the fast subsidence mode; (B) sketches showing the combinations of the pitching and horizontal motions in this mode. δu+q+ and δθ, disturbance quantities in non-dimensional x-component of velocity, pitching rate and pitch angle, respectively.

Fig. 9.

(A) Characteristic transients of disturbance quantities in the fast subsidence mode; (B) sketches showing the combinations of the pitching and horizontal motions in this mode. δu+q+ and δθ, disturbance quantities in non-dimensional x-component of velocity, pitching rate and pitch angle, respectively.

### The slow subsidence mode

For this mode, thalf is 57.8(Table 4); it takes about 58 wingbeats for the disturbance to decrease to half of its initial value. Unlike the above two modes, in which δw+ is negligibly small, this mode is a motion in which δw+ is the main variable (Table 5); other variables are one order of magnitudes or more smaller. Since

$${\delta}{\dot{z}}_{\mathrm{E}}^{+}={\delta}w^{+}$$
(Equation 9), this mode represents a descending (or ascending) motion, with the descending (or ascending) rate decreasing relatively slowly [after 200 wingbeats, the descending (or ascending) rate decreases to 5% of its initial value]. The characteristic transient of δw+ is plotted in Fig. 10.

Fig. 10.

Characteristic transients of disturbance quantities in the slow subsidence mode; δw+, disturbance quantities in non-dimensional z-component of velocity.

Fig. 10.

Characteristic transients of disturbance quantities in the slow subsidence mode; δw+, disturbance quantities in non-dimensional z-component of velocity.

### Physical interpretation of the motions of the natural modes

The eigenvalue and eigenvector analysis of the longitudinal small disturbance equations of the bumblebee has identified one unstable oscillatory mode and two (stable) monotonic subsidence modes. It is desirable to examine the physical processes of the motions of the natural modes and interpret the motions physically.

#### The unstable oscillatory mode

Let us examine the first cycle of the motion. For clarity, the first half-cycle of the characteristic transients of δu+(

$${\delta}{\dot{x}}_{\mathrm{E}}^{+}$$
⁠),δ q+ and δθ is replotted in Fig. 11. The physical process of the motion is sketched in Fig. 12.

Fig. 11.

Characteristic transients of disturbance quantities (first half-cycle) in the unstable oscillatory mode. δu+q+ and δθ, disturbance quantities in non-dimensional x-component of velocity, pitching rate and pitch angle, respectively.

Fig. 11.

Characteristic transients of disturbance quantities (first half-cycle) in the unstable oscillatory mode. δu+q+ and δθ, disturbance quantities in non-dimensional x-component of velocity, pitching rate and pitch angle, respectively.

Fig. 12.

Diagram of the bumble attitude, horizontal speed and pitching rate, and the forces and moments during the first cycle of the disturbed motion (unstable oscillatory mode). δu+q+ and δθ, disturbance quantities in non-dimensional x-component of velocity, pitching rate and pitch angle, respectively; ΔM+, non-dimensional pitching moment; mg, insect weight; F0, total aerodynamic force at equilibrium.

Fig. 12.

Diagram of the bumble attitude, horizontal speed and pitching rate, and the forces and moments during the first cycle of the disturbed motion (unstable oscillatory mode). δu+q+ and δθ, disturbance quantities in non-dimensional x-component of velocity, pitching rate and pitch angle, respectively; ΔM+, non-dimensional pitching moment; mg, insect weight; F0, total aerodynamic force at equilibrium.

At the beginning of the cycle (t+=0; Fig. 11), δθ is at its local minimum value, δu+(

$${\delta}{\dot{x}}_{\mathrm{E}}^{+}$$
⁠)is positive, and δq+ is zero; that is, at the beginning the bumblebee has a negative δθ and is moving forward with zero pitching rate. As seen in Fig. 12A, the negative δθ tilts forward the resultant aerodynamic force of reference flight (denoted as F0). The horizontal component of F0 tends to accelerate the forward motion (increasing δu+). The forward motion in turn produces a nose-up pitching moment (denoted by ΔM+;
$${\Delta}M^{+}=M_{\mathrm{u}}^{+}{\delta}u^{+}$$
and
$$M_{\mathrm{u}}^{+}$$
is positive), which would produce a nose-up pitching rate (δq+), making the magnitude of δθ to decrease (δθ to increase).

When δθ has increased to zero at t+≈10(the bumblebee has moved to the configuration shown in Fig. 12B),δ q+ does not reach its local maximum value, but continues to increase (see Fig. 11). This is because at this time δu+ is still large and so is the nose-up pitching moment(ΔM+). As a result, δθ would increase with time at a faster rate than when it is smaller than zero, which would cause the amplitude of δθ to become larger than that in the preceding quarter cycle. We thus see that the combination of forward motion and nose-up pitching causes the instability.

Now δθ has become positive and F0 is tilted backwards, which would slow the forward motion. At t+≈18, δu+ decreases to zero and changes sign (the bumblebee has moved to the configuration of Fig. 12C). Then, the bumblebee moves backward. The backward motion would produce a nose-down pitching moment,reducing the nose-up pitching rate (δq+). At t+≈24, δq+ changes sign,δθ reaches its local maximum value (the bumblebee has moved to the configuration of Fig. 12D).

In the next half-cycle, the above process repeats in an opposite direction(Fig. 12D-A); here it is the combination of backward motion and nose-down pitching that produces the destabilizing effect.

#### The fast subsidence mode

In this mode, as seen in Table 5 and Fig. 9A, whenδθ and δu+ have positive initial values,δ q+ has a negative initial value. The positiveδθ tilts F0 backwards, the horizontal component of which tends to reduce δu+; the forward motion (δu+) produces a nose-up pitching moment(

$${\Delta}M^{+}{\approx}M_{\mathrm{u}}^{+}{\delta}u^{+}$$
⁠)that tends to reduce the nose-down pitching rate(δq+); the nose-down pitching rate(δq+) tends to reduce δθ. This results in the monotonic decay of the disturbance quantities. (When δθ andδ u+ have negative initial values,δ q+ has a positive one; the motion can be explained similarly.)

#### The slow subsidence mode

In this mode, when the bumblebee descends initially due to some disturbance(i.e. δw+ has a positive initial value), the descending rate will decrease with time(Fig. 10). A positiveδ w+ produce a upward force(

$${\Delta}Z^{+}{\approx}Z_{\mathrm{w}}^{+}{\delta}w^{+}$$
⁠;
$$Z_{\mathrm{w}}^{+}$$
is negative), which tends to decrease the descending rate δw+, and the decreasedδ w+ would in turn reduce the upward force, resulting in the monotonic decay of the descending rate. (Whenδ w+ has a negative initial value, the corresponding motion can be explained similarly.)

#### The effects of the rate derivatives

In the preceding discussion, the effects of the rate derivatives(

$$X_{\mathrm{q}}^{+}$$
⁠,
$$Z_{\mathrm{q}}^{+}$$
and
$$M_{\mathrm{q}}^{+}$$
) are not mentioned. Here we discuss their effects on the motions. Since the reference flight is hovering flight,the magnitudes of
$$X_{\mathrm{q}}^{+}$$
and
$$Z_{\mathrm{q}}^{+}$$
are close to zero and
$$M_{\mathrm{q}}^{+}$$
(=-0.88) is relatively large (see Table 1). It is expected that
$$M_{\mathrm{q}}^{+}$$
has provided damping effect to the system; that is, without
$$M_{\mathrm{q}}^{+}$$
, the unstable oscillatory mode would grow faster and the stable subsidence modes would decay more slowly. In order to see this quantitatively, we set the rate derivatives in A to zero and computed the corresponding eigenvalues and eigenvetors. The results are shown in Tables 6 and 7, respectively. Comparing the results in Tables 6 and 7 with those in Tables 2 and 5, we see that without the damping effect, the growth rate of the oscillatory mode is 64% larger and the decaying rate of the fast subsidence mode is 20% smaller than those in the case with the damping effect. It is interesting to note that the growth rate of the slow subsidence mode is the same with or without the damping effect. This is because in this mode, δq+ is negligibly small.

Table 6.

Eigenvalues of the system matrix when the rate derivatives are set to zero

Mode 1, λ1,2Mode 2, λ3Mode 3, λ4
0.074±0.134i -0.158 -0.012
Mode 1, λ1,2Mode 2, λ3Mode 3, λ4
0.074±0.134i -0.158 -0.012

λ1,2, a pair of complex conjugate eigenvalues.λ 3 and λ4, real eigenvalues.

i, imaginary number.

Table 7.

Magnitudes and phase angles of the components of each of the three eigenvectors when the rate derivatives are set to zero

Modeδu+δw+δq+δθ
Unstable oscillatory 0.9×10-1 (121.7°) 1.9×10-4 (64.6°) 1.5×10-1 (60.9°) 1 (0°)
Fast subsidence 0.9×10-1 (0°) 2.2×10-4 (180°) 1.6×10-1 (180°) 1 (0°)
Slow subsidence 1.4 (0°) 1.7×10 (0°) 1.2×10-2 (180°) 1 (0°)
Modeδu+δw+δq+δθ
Unstable oscillatory 0.9×10-1 (121.7°) 1.9×10-4 (64.6°) 1.5×10-1 (60.9°) 1 (0°)
Fast subsidence 0.9×10-1 (0°) 2.2×10-4 (180°) 1.6×10-1 (180°) 1 (0°)
Slow subsidence 1.4 (0°) 1.7×10 (0°) 1.2×10-2 (180°) 1 (0°)

δu+, δw+q+ and δθ, disturbance quantities in non-dimensional x-component and z-component of velocity,pitching rate and pitch angle, respectively.

Numbers in the parentheses are the phase angles.

### The rigid body approximation

Taylor and Thomas (2002)have discussed the constraints on the rigid body approximation in detail. In general, the rigid body approximation only works well if the wingbeat frequency is at least an order of magnitude (10 times) higher than the highest frequency of the natural modes. They reasoned, using reduced order approximations to the natural modes of motion, that this could be expected to be true in animal flight. In the present study on the disturbed longitudinal motion of the hovering bumblebee, the period of the oscillatory mode is about 50 times the wingbeat period (see Table 4), which is much more than 10 times larger than the wingbeat period.

It should be noted that in the fast subsidence mode, a disturbance quantity varies from its initial (maximum) value to half of the value in 3.5 wingbeats. Is this too short to apply the rigid body approximation? To answer this question, let us look at an oscillatory mode, in which a variable would usually vary from a peak value to half of the value in a time 16% of its period. For the rigid body approximation to be appropriate, as stated above,the period should be at least 10 times as long as the wingbeat period; then 16% of the period is at least 1.6 wingbeat period. We thus see that in the fast subsidence mode, a variable decreasing from its initial value to half of that value in 3.5 wingbeat periods is slow enough for the rigid body approximation.

The above discussion shows that application of the rigid body approximation in the present analysis is appropriate. The result here provides an example that supports the reasoning of Taylor and Thomas(2002) on the applicability of the rigid body approximation to animal flight.

### The inherent dynamic stability and the equilibrium flight

The present model simulates the inherent dynamic stability of the bumblebee in the absence of active control. That is, in the disturbed motion, the model bumblebee uses the same wing kinematics as in the reference flight. A real bumblebee, if the motion is dynamically stable and the disturbances die out fast, might not make any adjustment to its wing kinematics, and could return to the equilibrium automatically'. In this case, the disturbed motion history predicted by the model represents that of the real bumblebee. In general,however, a real bumblebee makes continuous adjustments to its wing kinematics in order to keep to the reference flight, and the disturbed motion predicted by the model would be altered at an early stage.

In the present study, some of wing kinematic parameters (n, Φ,etc.) at reference flight are taken from or determined from the experimental data of Dudley and Ellington(1990a) and Ellington(1984b), and the others(αu, αd and

$${\bar{{\phi}}}$$
⁠) are solved from the force and moment equilibrium conditions. Because some simplifications are made in the model (e.g. the wing is a rigid flat plate, the translational velocity of the wing varies according to the simple harmonic function, etc.), the kinematic parameters that have been solved are only an approximation to those actually used by the bumblebee. Therefore, equilibrium flight, the dynamic stability of which our model studies, is only an approximation to the actual equilibrium flight of the bumblebee. Here, we must assume that the stability properties obtained by the model therefore only apply to the actual equilibrium flight of the bumblebee.

### The flight is unstable but the growth of the disturbances is relatively slow

As mentioned above, in general the disturbed motion is a linear superposition of the simple motions represented by the natural modes. For the hovering bumblebee, when disturbed from its reference flight, the disturbed motion is a linear combination of an unstable oscillatory mode and two stable subsidence modes. The growth of the disturbed motion is determined by the unsteady oscillatory mode. The unstable oscillatory mode doubles its amplitude in about 15 wingbeats (Table 4), which is about 0.1 s (the wingbeat period is 6.4 ms). To a person or a man-made machine, this growth rate is fast. But to a bumblebee,which can change its wing motions within a fraction of a wingbeat period, this growth rate might not be fast; the insect, if wishing to keep to the reference flight, has plenty of time to adjust its wing motion before the disturbances have grown large. For example, in the forward moving phase(Fig. 12A-C), the bumblebee might slightly decrease and increase the angle of attack of the wings during the downstroke and during the upstroke, respectively, from the equilibrium value of the angle of attack, and in the backward moving phase(Fig. 12D-F), the bumblebee might do the opposite. This would produce effects on the lift and the drag opposite to those produced by δu+, thus the destabilizing δM+ could be eliminated.

List of symbols

• A

system matrix

•
• c

mean chord length

•
• CL,w

vertical force coefficient of wing

•
• L,w

mean vertical force coefficient of wing

•
• CM,w

pitching moment coefficient of wing

•
• M,w

mean pitching moment coefficient of wing

•
• CT,w

thrust coefficient of wing

•
• T,w

mean thrust coefficient of wing

•
• Db

body drag

•
• e

reference flight condition

•
• E

earth

•
• F0

aerodynamic force of reference flight

•
• g

the gravitational acceleration

•
• i

imaginary number,

$$i=\sqrt{-1}$$

•
• Ib

pitching moment of inertia of the body about wing-root axis

•
• Iy

pitching moment of inertia about the y-axis of insect body

•
• l

length

•
• λ

distance from anterior tip of body to center of mass divided by body length

•
• lb

body length

•
• λl

distance from wing base axis to center of mass divided by body length

•
• λb

body length divided by R

•
• Lb

body lift

•
• Lw

vertical force of wing

•
• m

mass of the insect

•
• M

total aerodynamic pitching moment about center of mass

•
• M+

non-dimensional total aerodynamic pitching moment about center of mass

•
• $$M_{\mathrm{q}}^{+}$$

derivative of M+ with respect to q+

•
• $$M_{\mathrm{u}}^{+}$$

derivative of M+ with respect to u+

•
• $$M_{\mathrm{w}}^{+}$$

derivative of M+ with respect to w+

•
• n

stroke frequency

•
• o, o′, o1, oE

origins of the frames of reference

•
• q

pitching angular-velocity about the center of mass

•
• q+

non-dimensional pitching angular-velocity about the center of mass

•
• r2

radius of the second moment of wing area

•
• R

wing length

•
• Re

Reynolds number

•
• S

area of one wing

•
• St

area of two wings

•
• t

time

•
• tdouble

time for a divergent motion to double in amplitude

•
• thalf

time for a divergent motion to half in amplitude

•
• tw

period of the wingbeat cycle

•
• non-dimensional time (=0 at the start of a downstroke and =1 at the end of the subsequent upstroke)

•
• T

period of the insect motion

•
• Tw

horizontal force of wing

•
• u

component of velocity along x-axis

•
• u+

component of non-dimensional velocity along x-axis

•
• w

component of velocity along z-axis

•
• w+

component of non-dimensional velocity along z-axis

•
• ut

translational velocity of the wing

•
• ut+

non-dimensional translation velocity of the wing

•
• U

reference velocity

•
• x, y, z

coordinates in the body-fixed frame of reference (with origin at center of mass)

•
• x′, y′, z′

coordinates in the frame of reference with origin at wing root and z′ perpendicular to stroke plane

•
• x1, y1, z1

coordinates in the frame of reference with origin at wing root and z1 in vertical direction

•
• xE, yE, zE

coordinates in a system fixed on the earth

•
• E

xE-component of the velocity of the mass center of the insect

•
• X

x-component of the total aerodynamic force

•
• X+

non-dimensional x-component of the total aerodynamic force

•
• $$X_{\mathrm{q}}^{+}$$

derivative of X+ with respect to q+

•
• $$X_{\mathrm{u}}^{+}$$

derivative of X+ with respect to u+

•
• $$X_{\mathrm{w}}^{+}$$

derivative of X+ with respect to w+

•
• E

zE-component of the velocity of the mass center of the insect

•
• Z

z-component of the total aerodynamic force

•
• Z+

non-dimensional z-component of the total aerodynamic force

•
• $$Z_{\mathrm{q}}^{+}$$

derivative of Z+ with respect to q+

•
• $$Z_{\mathrm{u}}^{+}$$

derivative of Z+ with respect to u+

•
• $$Z_{\mathrm{w}}^{+}$$

derivative of Z+ with respect to w+

•
• α

geometric angle of attack of wing

•
• $${\dot{{\alpha}}}$$

angular velocity of pitching rotation

•
• $${\dot{{\alpha}}}^{+}$$

non-dimensional angular velocity of pitching rotation

•
• $${\dot{{\alpha}}}_{0}^{+}$$

a constant

•
• αd

midstroke geometric angle of attack in downstroke

•
• αu

midstroke geometric angle of attack in upstroke

•
• β

stroke plane angle

•
• δ

small disturbance notation (prefixed to a perturbed state variable)

•
• Δ

increment notation

•
• θ

pitch angle between the x-axis and the horizontal

•
• λ

generic notation for an eigenvalue

•
• Δλr

duration of wing rotation or flip duration (non-dimensional)

•
• ρ

density of fluid

•
• τ

non-dimensional time

•
• τc

non-dimensional period of one flapping cycle

•
• τr

non-dimensional time when pitching rotation starts

•
• ν

kinematic viscosity

•
• φ

azimuthal or positional angle

•
• $${\bar{{\phi}}}$$

mean positional angle

•
• $${\dot{{\phi}}}$$

angular velocity of azimuthal rotation

•
• $${\dot{{\phi}}}^{+}$$

non-dimensional angular velocity of azimuthal rotation

•
• Φ

stroke amplitude

•
• χ

body angle

•
• χ0

free body angle

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

Dickinson, M. H. and Götz, K. G. (
1993
). Unsteady aerodynamic performance of model wings at low Reynolds numbers.
J. Exp. Biol.
174
,
45
-64.
Dickinson, M. H., Lehman, F. O. and Sane, S. P.(
1999
). Wing rotation and the aerodynamic basis of insect flight.
Science
284
,
1954
-1960.
Dudley, R. and Ellington, C. P. (
1990a
). Mechanics of forward flight in bumblebees. I. Kinematics and morphology.
J. Exp. Biol.
148
,
19
-52.
Dudley, R. and Ellington, C. P. (
1990b
). Mechanics of forward flight in bumblebees. II. Quasi-steady lift and power requirements.
J. Exp. Biol.
148
,
53
-88.
Ellington, C. P. (
1984a
). The aerodynamics of hovering insect flight. II. Morphological parameters.
Phil. Trans. R. Soc. Lond. B
305
,
17
-40.
Ellington, C. P. (
1984b
). The aerodynamics of hovering insect flight. III. Kinematics.
Phil. Trans. R. Soc. Lond. B
305
,
79
-113.
Ellington, C. P., van den Berg, C. and Willmott, A. P.(
1996
). Leading edge vortices in insect flight.
Nature
347
,
472
-473.
Etkin, B. (
1972
).
Dynamics of Atmospheric Flight
. New York: John Wiley and Sons,Inc.
Sun, M. and Tang, J. (
2002
). Unsteady aerodynamic force generation by a model fruit fly wing in flapping motion.
J. Exp. Biol.
205
,
55
-70.
Sun, M. and Wu, J. H. (
2003
). Aerodynamic force generation and power requirements in forward flight in a fruit fly with modeled wing motion.
J. Exp. Biol.
206
,
3065
-3083.
Taylor, G. K. and Thomas, A. L. R. (
2002
). Animal flight dynamics. II. Longitudinal stability in flapping flight.
J. Theor. Biol.
214
,
351
-370.
Taylor, G. K. and Thomas, A. L. R. (
2003
). Dynamic flight stability in the desert locust Schistocerca gregaria.
J. Exp. Biol.
206
,
2803
-2829.
Thomas, A. L. R. and Taylor, G. K. (
2001
). Animal flight dynamics. I. Stability in gliding fight.
J. Theor. Biol.
212
,
399
-424.
Usherwood, J. R. and Ellington, C. P. (
2002a
). The aerodynamics of revolving wings. I. Model hawkmoth wings.
J. Exp. Biol.
205
,
1547
-1564.
Usherwood, J. R. and Ellington, C. P. (
2002b
). The aerodynamics of revolving wings. II. Propeller force coefficients from mayfly to quail.
J. Exp. Biol.
205
,
1565
-1576.
Wang, Z. J. (
2000
). Two dimensional mechanism for insect hovering.
Phys. Rev. Lett.
85
,
2216
-2219.
Wu, J. H. and Sun, M. (
2004
). Unsteady aerodynamic forces of a flapping wing.
J. Exp. Biol.
207
,
1137
-1150.