First published online January 25, 2005
Journal of Experimental Biology 208, 447-459 (2005)
Published by The Company of Biologists 2005
doi: 10.1242/jeb.01407
Dynamic flight stability of a hovering bumblebee
Mao Sun* and
Yan Xiong
Institute of Fluid Mechanics, Beijing University of Aeronautics and
Astronautics, Beijing 100083, People's Republic of China
*
Author for correspondence (e-mail:
m.sun{at}263.net)
Accepted 24 November 2004
 |
Summary
|
|---|
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.
Key words: dynamic stability, flapping flight, hovering, bumblebee, insect, Navier-Stokes simulation, natural modes of motion, bumblebee
|
Introduction
|
|---|
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 Gö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.
|
Materials and methods
|
|---|
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.
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 (
).
oExEyEzE
is 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:
 | (1) |
 | (2) |
 | (3) |
 | (4) |
 | (5) |
 | (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+Xuu+Xww+Xqq,
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:
 | (7) |
 | (8) |
 | (9) |
where A is the system matrix:
 | (10) |
where
am=0.5
USttw,
and ag = U/tw, and the non-dimensional
forms are: u+ = u/U,
w+ = w/U, q+ =
qtw, X =
X/0.5U2St ( denotes
the air density and St denotes the area of two
wings),Z+=Z/0.5U2St,
M+=Z/0.5U2Stc,
t+=t/tw,
and
(using the flight data given below, am,
aI and ag are computed as
am=1.96 mg,
aI=0.233x10-9 kg m2,
ag=710.9 m s-1; is 1.25 kg m-3
and 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.48x10-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
. 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).
View larger version (96K):
[in this window]
[in a new window]
|
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
(
) is related to
,
where 
is non-dimensional time. For the bumblebee, on the basis of data
given by Dudley and Ellington
(1990a), ut
is approximated by the simple harmonic function:
 | (11) |
where the non-dimensional translational velocity
,
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 d
for the downstroke and u for the upstroke. Around stroke
reversal, the wing flips and changes with time. The angular velocity
(
) is given by:
 | (12) |
where the non-dimensional form
,
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 d,
u and r are specified,
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.
View larger version (23K):
[in this window]
[in a new window]
|
Fig. 3. Sketches of the reference frames and wing motion.
o1x1y1z1
is an inertial frame, with the x1,y1
plane in the horizontal plane.
o'x'y'z' is
another inertial frame, with the x'y' plane in
the stroke plane.  , positional angle of the wing; min
and 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
(
) 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.22c and the flip
rotation is assumed to be symmetrical (thus the flip timing r
is determined in terms of r). d,
u and 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 Tw
not 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:
 | (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.5U2St
and 0.5U2Stc,
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
(). In the present study,
d, u and
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.5U2St, 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 lxl real
matrix has l eigenvalues (
1,
2,... 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:
 | (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:
 | (15) |
A pair of complex conjugate eigenvalues, e.g.
,
will result in oscillatory time variation of the disturbance quantities with
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:
 | (16) |
and the times to double or half the oscillatory amplitude are
 | (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.
|
Results
|
|---|
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 93x109x78 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 71x73x96
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
, the mean vertical and horizontal
forces and mean pitching moment of the wings would be different.
d, u and
are determined using equilibrium
conditions. The calculation proceeds as follows. A set of values for
d, u and
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
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
=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
,
,
,
,
,
,
,
and
, estimated using the data in
Fig. 4, are shown in
Table 1.
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, is almost zero,
is negative and
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,w,
CT,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
. 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
. 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 .
View larger version (21K):
[in this window]
[in a new window]
|
Fig. 5. Time courses of CL,w (A),
CT,w (B) and CM,w (C) in
one flapping cycle. CL,w,
CT,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+.
and
are very small and
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 ;
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 and
. This explains the relatively large
and small
and
.
View larger version (19K):
[in this window]
[in a new window]
|
Fig. 6. Time courses of CL,w (A),
CT,w (B) and CM,w (C) in
one flapping cycle. CL,w,
CT,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,
and
are very small and
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 and
. 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
. 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.
View larger version (17K):
[in this window]
[in a new window]
|
Fig. 7. Time courses of CL,w (A),
CT,w (B) and CM,w (C) in
one flapping cycle. CL,w,
CT,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.
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.
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.
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
(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+
(
),
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.
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.
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
(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.
View larger version (6K):
[in this window]
[in a new window]
|
Fig. 10. Characteristic transients of disturbance quantities in the slow subsidence
mode; w+, disturbance quantities in non-dimensional
z-component of velocity.
|
|
|
Discussion
|
|---|
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+
(),
q+ and is replotted in
Fig. 11. The physical process
of the motion is sketched in Fig.
12.
At the beginning of the cycle (t+=0;
Fig. 11), is at
its local minimum value, u+
()
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+;
and 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
(
)
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
(
;
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
(,
and
) are not mentioned. Here we discuss
their effects on the motions. Since the reference flight is hovering flight,
the magnitudes of and
are close to zero and
(=-0.88) is relatively large (see
Table 1). It is expected that
has provided damping effect to the
system; that is, without , 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.
TOP
Summary
Introduction
