## SUMMARY

We combine three-dimensional descriptions of the movement patterns of the
shoulder, elbow, carpus, third metacarpophalangeal joint and wingtip with a
constant-circulation estimation of aerodynamic force to model the wing
mechanics of the grey-headed flying fox (*Pteropus poliocephalus*) in
level flight. Once rigorously validated, this computer model can be used to
study diverse aspects of flight. In the model, we partitioned the wing into a
series of chordwise segments and calculated the magnitude of segmental
aerodynamic forces assuming an elliptical, spanwise distribution of
circulation at the middle of the downstroke. The lift component of the
aerodynamic force is typically an order of magnitude greater than the thrust
component. The largest source of drag is induced drag, which is approximately
an order of magnitude greater than body form and skin friction drag. Using
this model and standard engineering beam theory, we calculate internal
reaction forces, moments and stresses at the humeral and radial midshaft
during flight. To assess the validity of our model, we compare the
model-derived stresses with our previous *in vivo* empirical
measurements of bone strain from *P. poliocephalus* in free flapping
flight. Agreement between bone stresses from the simulation and those
calculated from empirical strain measurements is excellent and suggests that
the computer model captures a significant portion of the mechanics and
aerodynamics of flight in this species.

## Introduction

Although the link between organismal structure and function has been recognized for centuries, it is clear that the strength of this relationship varies tremendously. Evolutionary theory suggests that this relationship should be particularly strong when the consequences of deviations from optimal designs are energetically costly. In this light, the relationship between the structural design of the wing and the mechanics and energetics of flight in the Order Chiroptera constitutes an informative case study for several reasons: (i) powered flight imposes extreme mechanical and energetic demands on the locomotor system (Kurta et al., 1989⇓; Swartz, 1997⇓; Swartz, 1998⇓; Thomas, 1975⇓; Thomas and Suthers, 1972⇓; Winter et al., 1993⇓), raising the metabolic cost of deviations from optimal mechanical design relative to that for terrestrial locomotion in quadrupeds; (ii) dermopterans and primates, the nearest non-flying relatives of bats (Simmons, 1995⇓), provide considerable comparative information about the likely nature of bat ancestors; and (iii) the diversity in body size, flight style, wing shape and phylogenetic affinity of the nearly 1000 extant species of bats furnish rich comparative material within the group.

Both classic and recent results demonstrate that the bat wing is unique among mammalian limbs in anatomical design and mechanical function and suggest that the specialized features of the bat musculoskeletal system are linked directly to flight capabilities (Findley et al., 1972⇓; Hermanson and Alternbach, 1983⇓; Holbrook and Odland, 1978⇓; Norberg, 1970a⇓; Norberg, 1972a⇓; Norberg, 1972b⇓; Papadimitriou et al., 1996⇓; Strickler, 1978⇓; Swartz, 1997⇓; Swartz, 1998⇓; Swartz et al., 1992⇓; Swartz et al., 1996⇓; Vaughan, 1959⇓; Vaughan, 1970b⇓). Similarly, significant contributions have been made to our understanding of the kinematics of bat flight (e.g. Aldridge, 1986⇓; Aldridge, 1987⇓; Baggøe, 1987⇓; Brandon, 1979⇓; Norberg, 1970b⇓; Norberg, 1976a⇓; Norberg 1976b⇓; Norberg, 1990⇓; Rayner, 1987⇓; Rayner et al., 1986⇓; Vaughan, 1970a⇓). Some components of these kinematic and morphological studies have also contributed directly to our understanding of flight mechanics, including the interrelationship between wing membrane tension and aerodynamic force (Norberg, 1972a⇓; Pennycuick, 1973⇓), the velocity-dependence of wing kinematics (Aldridge, 1986⇓) and the relationship between wing inertia, energetics and maneuverability (Norberg, 1976a⇓; Thollesson and Norberg, 1991⇓).

However, despite decades of study, we understand much less about the mechanics and energetics of flying than we do of walking, running or swimming in vertebrates; this, in turn, limits our knowledge of how flight capacity and the intricate array of physiological and morphological specializations associated with it have evolved. In part, this is because our insight is constrained by our limited understanding of the nature of the complex and dynamically changing forces experienced by the wing during flapping flight and the technical difficulties of approaching this subject empirically. To meet these challenges and thereby gain new insight into the mechanics and energetics of the bat wing, we have developed a three-dimensional computer simulation of bat flight. Computer modeling is a powerful approach to gaining insight into the complexities of animal flight and has been adopted with success in diverse ways by students of bat, bird and insect flight (e.g. Bennett, 1977⇓; DeLaurier, 1993a⇓; Ellington, 1975⇓; Ellington, 1978⇓; Ellington, 1984⇓; Ennos, 1988⇓; Norberg, 1975⇓; Norberg, 1970b⇓; Norberg, 1976a⇓; Pennycuick, 1968⇓; Pennycuick, 1975⇓; Rayner, 1986⇓; Rees, 1975⇓; Spedding, 1992⇓; Spedding and DeLaurier, 1995⇓; Withers, 1981⇓).

Our model computes wing bone stresses, joint forces and moments and other mechanical and energetic parameters from wing kinematics and structural geometry placed in a context of a well-founded, realistic and detailed aerodynamic model. The wing anatomy of bats is particularly well suited to this kind of engineering analysis: the wing comprises a jointed network of virtually rigid structural supports interconnected by an essentially two-dimensional elastic membrane; bats operate at Reynolds numbers high enough for appropriate application of inviscid aerodynamic theory; and bat flight occurs at Strouhal numbers at which complex unsteady forces are far less important than for insect flight. In this context, our model is structured to estimate accurately the forces experienced by the wing, to analyze in detail one of these forces, the internal forces developed in the wing, and, ultimately, to detail further the contributions to a critical element of the internal forces, the wing membrane forces.

In addition to its detail and accuracy, this model is noteworthy in that we
are able to validate key aspects by direct comparison of the model's estimates
of wing bone stresses with empirically measured values (Swartz et al.,
1992⇓). A well-validated model
of bat flight should be able to reproduce, in order of increasing
sophistication, appropriate skeletal loading (tension *versus*
compression), the general pattern of change in skeletal stress in relation to
the wingbeat cycle, stresses of realistic magnitude and details of changes in
stresses during the wingbeat. We are also able to compare the predictions of
the model concerning the vertical movements of the animal's center of mass
with measurements made directly from wind-tunnel flights (Bartholomew and
Carpenter, 1973⇓).

Here, we describe the structure of the model, evaluate its ability to reproduce important aspects of bat flight mechanics realistically and examine the sensitivity of the model to various inputs and assumptions. Once validated and evaluated for parametric sensitivity, we intend to apply the model to a broader comparative analysis of the flight mechanics of bats that differ in wing morphology and/or flight behavior. We believe this model can be employed to gain insight into diverse problems in the mechanics and evolution of bat flight in future studies.

## Materials and methods

### Model outline

Our model of bat flight is based on the grey-headed flying fox,
*Pteropus poliocephalus* Temminck. The model is constructed to employ
empirical descriptions of morphology and kinematics, along with reasonable
assumptions concerning aerodynamics, to estimate accurately all but one of the
forces experienced by a bat wing (see below). The model then solves for the
remaining force, the internal force carried by the wing structures. For this
analysis, we subdivided the wing into a series of segments or strips and
measured the shape and mass of each wing segment
(Fig. 1). Instantaneous force
balance was then applied locally to each wing segment or strip. We quantified
the three-dimensional motion of critical wing landmarks over the wingbeat
cycle using high-speed films of a *P. poliocephalus* made in a wind
tunnel (Bartholomew and Carpenter,
1973⇓) and used this information
to compute the inertial forces of each wing segment (see below, section on
dynamics). We modeled the force induced by air flow around the wing as an
aerodynamic force, separated into lift and drag components, plus an added mass
force. In this model of level flight at constant velocity, we assumed constant
circulation (see below, section on aerodynamic force distribution and
constant-circulation flight). Because of this assumption, it is unnecessary to
quantify wing segment camber, incident air speed, angle of attack and unsteady
aerodynamic effects. We also assumed an elliptical, spanwise distribution of
aerodynamic force (see below, section on aerodynamic force distribution and
constant-circulation flight) and adjusted the amplitude of the aerodynamic
force to balance mean lift with body weight over one wingbeat cycle.

We employed Newton's second law of motion to balance the inertial force of a wing strip with the external forces as well as the internal force within the wing segment, and then used the segmental balance of forces at each instant in time to solve for the internal force carried by the wing structures. We calculated the stresses developed at the midshaft of the humerus and radius by summing the internal forces over all strips distal to the bone site of interest and computing the moments of these forces about the bone's midshaft. We modified these maximum possible internal forces and moments by subtracting the forces and moments transmitted directly by the plagiopatagium (armwing) to the body, and then employed standard engineering beam theory to convert the remaining three orthogonal forces and moments at the midshaft into normal and shear stresses at the bone surface, and compared stresses calculated for the humeral and radial midshafts with stresses inferred from strains measured empirically in flying bats (Swartz et al., 1992⇓). The close correspondence between the model estimates and the empirical data provides robust evidence that the model approximates the mechanics of bat flight in biologically meaningful ways. As further validation, we also compared the vertical displacement of the animal's center of mass during the wingbeat cycle, as predicted by the model, with the empirically observed vertical motion.

The computer program that solves the equations and algorithms presented
below was written in Fortran 77 (Pro Fortran, version 5.0, Absoft) and run on
a Power Macintosh G3. The program provides outputs (bone stresses, joint
moments, etc.) at 40 time steps evenly spaced over the course of a single
wingbeat cycle. Hence, the time *t* corresponding to the integer value
of the time step *q* is given by:

since 40 time steps are used in the period of one wingbeat. Because the force balance employed in the model applies at every instant of flight, the discretized time step does not incur any approximation other than the possibility of interpolating model output between time steps.

### Choice of species

We selected the species *Pteropus poliocephalus*, the greyheaded
flying fox, for the first application of our bat flight model for several
reasons. First, this species occurs at relatively high densities in the
proximity of Brisbane, Australia, and the University of Queensland and has
therefore been a subject of previous research (e.g. Carpenter,
1985⇓; Swartz,
1998⇓; Swartz et al.,
1992⇓; Swartz et al.,
1993⇓). Second, individuals are
large (adult body mass typically 550-950 g), facilitating mechanical
assessment of functionally important wing structures, including *in
vivo* measurements of wing bone stresses. Therefore, we possess detailed
information concerning the *in vivo* loading of the humerus, radius,
metacarpals III and V and proximal phalanges III and V for this species.
Third, high-speed dorsal, lateral and oblique films of individuals of this
species flying in a wind tunnel provide detailed information concerning wing
kinematics (Bartholomew and Carpenter,
1973⇓). Fourth, level flight is
probably an ecologically relevant flight mode for this species given that
these fruit-eating bats often fly 30-50 km during foraging bouts on a single
night and may migrate hundreds of kilometers in a year (Eby,
1991⇓).

### Validation by in vivo strain recordings

Measurements of bone surface principal strain magnitudes and orientations
were made from eight individual *P. poliocephalus* in previous studies
(Swartz et al., 1992⇓; Swartz
et al., 1993⇓). In these
studies, animals were wild-caught and trained to fly the length of a 30 m
flight cage without stopping. Within 2 weeks of capture, rosette strain gauges
were surgically implanted on the subperiosteal surfaces of the midshafts of
wing bones, and the animals recovered fully from the effects of surgery. We
then collected data from up to nine strain gauge elements simultaneously
*via* a lightweight cable (100 Hz) and synchronized the data with video
recordings of the animals' wing movements. Data from individual rosette
elements were analyzed to obtain maximum and minimum principal strain
magnitudes and orientations, and strain values were converted to stresses
assuming that the compact cortical bone of the long bones of the wings of
large bats is similar in its mechanical properties to that of other mammals
and birds (Carter, 1978⇓; Beer
and Johnston, 1981⇓; Biewener,
1983⇓; Currey,
1987⇓). Model predictions of
bone stresses were then calculated for specific anatomical sites from which
strain data had been collected, and the experimentally determined stress
profiles for a given recording site were normalized to a standardized wingbeat
cycle, synchronizing the mid-downstroke, downstroke—upstroke transition
and upstroke—downstroke transition.

### Morphological parameters and general flight characteristics

Mass, wing size and shape and the dimensions of individual wing bones were
assessed by direct weight measurement, tracings of wing outlines and
measurements of high-resolution radiographs of a member of the study
population used in previous work (Swartz et al.,
1992⇓)
(Table 1). Although the model
data were measured from a single individual killed following completion of
bone strain recording, these parameters vary relatively little among
individuals of a given body mass, and the small intraspecific variation would
have no significant effect on model results. Video recordings of flying foxes
were used to estimate the flight speed of the bats (Swartz et al.,
1992⇓). A typical flight speed
*U* of approximately 6 m s^{-1} was estimated from the time
needed to fly approximately 25 m; this represents a moderate speed for this
species. We measured the wingbeat period and detailed kinematics (see also
section on wing kinematics) from a high-speed film of a flying fox in a wind
tunnel. Wing length was measured as the distance from the proximal shoulder to
the proximal carpus plus the distance from the carpus to the wingtip for a
fully extended handwing. Body width is taken as the shoulder-to-shoulder
distance. We calculated the mean area of a single wing from plan-view
photographs of a deeply anesthetized study subject with its wings fully
extended.

### Wing segmentation

In the computer model, we conceptually divided the wing at mid-downstroke
into a series of rectangular segments of variable chord and width, determined
anatomically (see below) (DeLaurier,
1993b⇓; Norberg,
1976a⇓; Thollesson and Norberg,
1991⇓). We subdivided the wing
into 14 rectangular segments: seven segments of equal width between the
shoulder and the carpus, three equal-width segments between the carpus and the
third metacarpophalangeal (MCP) joint and four equal-width segments between
the third MCP joint and the wingtip (Fig.
1). Because the wing changes its three-dimensional conformation
during the wingbeat, we calculated the instantaneous segment widths and
segment motions by linear interpolation of the distances between the carpus,
third MCP joint and wingtip positions, respectively. Because all segments
between two proximodistally adjacent landmarks are defined to have equal
widths, elbow flexion causes the width of the seven armwing segments to
decrease, although we assumed that the chord of each rectangular segment
remains constant. We write the width of a wing segment as
*w*_{pq}, where the subscript *p* denotes the wing
segment and increases from proximal to distal, and the subscript *q*
denotes the time step during the wingbeat cycle. Since the model does not
employ lift coefficients, the angle of attack of each wing segment is not
required as input to the model. As a result, we consider each wing segment to
be parallel to the cranial—caudal axis at all times. However, wing
segment pitching angle and angle of attack variations are modeled
intrinsically through the calculation of the spatial orientation of
aerodynamic forces over time.

We sectioned the wing and, for each strip, we measured mass,
*m*_{p}, chord along the strip's midline,
*c*_{p}, leading edge position (perpendicular distance from the
center of the leading edge to a reference line connecting the left and right
glenohumeral joints), *e*_{p}, and center of mass position with
respect to the shoulder-to-shoulder reference line, *d*_{p}
(e.g. Thollesson and Norberg,
1991⇓)
(Fig. 1,
Table 2). Values of
*e*_{p} are negative, indicating a position cranial to the
reference line. We determined the position of the center of mass of each
segment by attaching it to a rigid cardboard rectangle of known mass and
locating the center of mass of the cardboard—wing ensemble. We shaped
the wing to match the plan view of an individual wing in the middle of the
downstroke; at this point in the wingbeat cycle, the wing is nearly horizontal
and coplanar. For the wing kinematics employed in the model (see below), the
middle of the downstroke corresponds to the dimensionless time
*t/T*=0.38. During the middle of the downstroke, we determined the
midline leading edge positions of all wing segments relative to straight lines
adjoining wing landmarks in order to position wing segments in the
cranial—caudal direction. We calculated the position and motion of
coplanar wing segments by linear interpolation between adjacent wing
landmarks. We kept the distances from segmental leading edges to the straight
lines connecting adjacent wing landmarks constant over the entire wingbeat
cycle, in effect conserving the shape of the flapping wing. Wing segment
position determined where the forces acting on that wing segment were
applied.

### Wing kinematics

In this section, we outline our method of describing in three dimensions the motion of the carpus, the third MCP joint and the wingtip. The three-dimensional positions of these major wing landmarks over one wingbeat cycle suffice to capture the large-scale features of wing motion; for the purposes of our model, we define the beginning of downstroke as the onset of downward motion of the carpus.

### Origin and axes for movement description

The body of a flying bat accelerates and decelerates vertically during
`level' flight and cannot therefore be used as an inertial frame of reference
from which to measure the accelerations of the animal's wing. Here, instead,
we describe the motion of the wing relative to the mean position of the
glenohumeral joint averaged over the entire wingbeat cycle, which constitutes
the origin of an inertial coordinate system moving with the mean forward
flight speed of the bat. We fixed the origin of a second coordinate system to
the right glenohumeral joint and related the motion of this non-inertial
coordinate system to the inertial coordinate system
(Fig. 2). We defined our axes
such that the *x* axis points horizontally to the right of the
direction of flight for the right wing (distally along the wing at
mid-downstroke), the *y* axis points vertically upwards (in the dorsal
direction) and the *z* axis points in the opposite direction to flight
(in the caudal direction along the cranial—caudal axis). We denote the
absolute acceleration of the glenohumeral joint (or shoulder) as the global
acceleration **a** with respect to the non-inertial origin. The global
acceleration of the shoulder is not known *a priori*, but is instead
computed internally by the model through the force balance on the body. The
unit vectors giving the direction of each coordinate axis are **i** for the
*x* axis, **j** for the *y* axis and **k** for the
*z* axis. We model the left wing implicitly through symmetry across the
*yz*, or midsagittal, plane. This symmetry is a reasonable assumption
for forward flight and gliding, but would not be appropriate when modeling
more complex flight maneuvers.

Three-dimensional coordinates of wing landmarks were taken from high-speed
movies of wind-tunnel flight (Bartholomew and Carpenter,
1973⇓). We selected these films
for analysis of wing kinematics because of their far greater resolution than
the video recordings made during previous *in vivo* strain measurement
experiments (Swartz et al.,
1992⇓; Swartz et al.,
1993⇓). The films were
converted to video and analyzed with the Peak Performance Motion Analysis
System (Peak Performance Technologies, Englewood, CO, USA) to obtain detailed
information regarding the dynamically changing positions of the carpus, third
MCP joint and wingtip relative to the right shoulder. The spatial resolution
of the digitizing process corresponds to approximately ±2 cm, or
approximately ±5 % of the mean wing length and ±12 % of the mean
wing chord. From both head-on (*xy* plane) and lateral (*yz*
plane) views, approximately 170 sets of two-dimensional coordinates were
obtained for one complete wingbeat cycle. The right shoulder coordinates were
subtracted from each wing landmark position. We combined the two views into a
composite of wing motion by overlaying the *x* coordinates of wing
motion at the carpus, third MCP joint and wingtip. We illustrate
representative motion of the carpus and wingtip projected onto the *xy*
plane with respect to the shoulder (Fig.
3).

The wingbeat cycle of flying foxes, like that of some birds, can be
partitioned into approximately 65 % downstroke and 35 % upstroke. Sinusoidal
motion, necessarily 50 % downstroke and 50 % upstroke, is therefore not
appropriate to describe these wing motions, and one would have to resort to
Fourier series expansions to describe wing kinematics; we have instead
selected cyclic polynomials to describe wing motion realistically. We
curve-fitted the time change in the *x, y* and *z* positions of
the carpus, third MCP joint and wingtip relative to the right shoulder with
eighth-order polynomials using a least-squares algorithm (KaleidaGraph,
version 3.0.5, Abelbeck Software) (Table
3). For example, we wrote the position of the carpus along the
*x* axis as:

where
*m*_{0},..., *m*_{8} are the fitted
coefficients, and the wingbeat period *T* is written explicitly as part
of the polynomial coefficients. One wingbeat cycle lasted from
0<*t/T*≤1 in nondimensional time regardless of variations in the
wingbeat period. We ensured that the polynomial was periodic in
0<*t/T*≤1 by enforcing identity in position, velocity and
acceleration at *t*=0 and *t*=1. We chose to impose these three
continuity conditions by the following three relationships among the
coefficients:

leaving six degrees of
freedom *m*_{0},..., *m*_{5} with which to fit
the data. All *x* direction curve fits have correlation coefficients
*r*≥0.959 whereas all *y* direction curve fits have
correlation coefficients *r*≥0.993. Displacements in the *z*
direction were an order of magnitude smaller than those in the other two
directions and matched the accuracy of the position data in all three
directions, resulting in curve fits with correlation coefficients
*r*>0.44.

Wing landmark positions near *t/T*=1 were not precisely periodic
because of rounding errors in the curve-fit coefficients
(Table 3). Moreover, double
differentiation of position curves with continuous acceleration guarantees
continuity of the acceleration at end points but not continuity in the slope
of the acceleration. These are subtle effects that become very pronounced when
taking derivatives of the position data. We enforced periodicity in all
position curves by using the value at *t/T*=0 in place of the value
obtained from *t/T*=1. Continuity and smoothness in velocity and
acceleration profiles was obtained by applying Savitzky—Golay smoothing
filters to the position, velocity and acceleration curves (Press et al.,
1989⇓). The amount of smoothing
was carefully tested to have no perceptible effect on landmark positions and
yet still eliminate the growth of spurious discontinuities near end points. We
calculated wing landmark velocities from second-order finite differences of
the smoothed position curves and accelerations from second-order finite
differences of the smoothed velocity curves. Finally, we calculated the
position, velocity and acceleration of each wing segment from linear
interpolation between coplanar wing segments, as described in the previous
section.

### Dynamics

From the curve fits of position data, we calculated the magnitude, orientation and location of the gravitational, inertial, added mass and aerodynamic forces acting on each wing segment over the course of one wingbeat cycle. We then solved for the internal force carried by wing structures within a segment by invoking Newton's second law of motion at discrete instants of time. We summed the internal forces within all wing segments to calculate the instantaneous global acceleration of the bat shoulder. We corrected the magnitudes of the segmental aerodynamic forces in order to enforce the level-flight criterion (no net vertical acceleration of the body over an entire wingbeat cycle). Finally, we calculated reaction forces and moments at the humeral and radial midshafts by assuming that each bone carries the entire internal load of all more distal wing structures minus the load carried directly to the body by the plagiopatagium (see below, section on forces, moments and stresses on the wing skeleton).

### Gravitational force

The action of gravity on the wing of a large bat is typically an order of
magnitude smaller than the aerodynamic force encountered in level flight
because the former acts only on the wing mass while the latter supports the
entire body mass (e.g. Thollesson and Norberg,
1991⇓). Nevertheless, the
segmental gravitational force is given by:
6
where *m*_{p} is the segmental mass and **j** is a unit
vector of the global *y* axis pointing in the direction opposite to
gravitational acceleration, ** g**. The gravitational force acts at
the center of mass of each wing segment (see
Table 2 for the chordwise
positions of the centers of mass relative to the positions of the leading
edge). We assumed that the center of mass lies at midspan of each wing
segment.

Near the wingtip, the trapezoidal and triangular shape of the last few wing segments would cause the center of mass to be located more proximal to the shoulder than the midline of the wing segment. We did not include this effect because segment widths and masses near the wingtip are sufficiently small that the effect of any correction is negligible.

### Inertial force

Inertial forces on the wing vary with wing motion and thus can range in
magnitude from zero during episodes of gliding to values comparable with the
large aerodynamic forces exerted on the wing during flapping flight. By
definition, the inertial force also represents the net sum of all external
forces as well as the force carried by the wing structures for a given wing
segment. On the basis of segmental accelerations, the inertial force on wing
segment *p* at time *q* is given by:
7
where each component of a force is treated independently by the model. The
inertial force acts at the center of mass and in the same vector direction as
the acceleration of a wing segment. The horizontal acceleration
*s*_{sz.q} and the vertical acceleration
*s*_{sy.q} are the two components of the global acceleration at
each instant of time and must be solved for iteratively in the program. We
note that the Einstein summation convention is not implied in this paper by
the use of repeated subscripts.

### Added mass force

The added mass force resists the acceleration of the wing and is sometimes
considered part of the inertial force since it is proportional to the
acceleration perpendicular to the plane of a wing segment (DeLaurier,
1993b⇓). The magnitude of this
acceleration is a^{.}n and it is always aligned with the unit normal
vector **n**. We define the normal vector **n**_{pq} at any time
*q* to point away from the dorsal face of wing segment *p*. The
component *n*_{z.pq} of the normal vector **n**_{pq}
is always zero given that we assume, as a first approximation, that wing
segments remain parallel to the cranial—caudal axis of the body at all
times. We write the added mass force for wing segment *p* at time
*q* as:

where
*w*_{pq} is the mean width of wing segment *p* at time
*q*, ρ0 is the density of air, and the added mass coefficient,
*C*_{m}=0.9, is that of a thin plate of finite width entraining
the volume *V* of air contained within the cylinder whose width is that
of the wing segment and whose diameter is that of the segment chord,
*V*=(π*w*_{pq}*c*_{p}^{2})/4.
The added mass force opposes wing segment acceleration and acts through the
center of a wing segment at *w*_{pq}/2 and
*c*_{p}/2. The magnitude of the added mass force on a given
wing segment is comparable with that of the inertial force whenever
*m*_{p}≈(0.9π*w*_{pq}*c*_{p}^{2}ρ0)/4
or *m*_{p}≈0.6 g for *P. poliocephalus*.

### Aerodynamic force

The aerodynamic force consists of lift and thrust components that arise
from the rates at which air is drawn up ahead of the wing and pushed down
behind the wing. In our study species, the aerodynamic forces are the largest
and most significant forces experienced by the wing, although segmental
inertial forces during a flapping cycle can momentarily achieve the magnitude
of segmental aerodynamic forces (see Results). Together, the aerodynamic and
inertial forces account for most of the external load exerted on wing
structures. Skin friction along the surface of the wings, form drag on the
body and induced drag from vortices in the wake provide the total drag
experienced during flight. We combined the total drag and the animal's weight
to find the total aerodynamic force and distributed this force over all the
wing segments. In the absence of detailed wing profiles in flight, we
prescribed the aerodynamic force generated by each wing segment through an
elliptical, spanwise distribution of circulation (Norberg,
1990⇓) at the middle of the
downstroke as a plausible approximation for a bat wing in fast forward flight
(Rayner et al., 1986⇓; Rayner,
1986⇓; Spedding,
1987⇓) (see also section on
aerodynamic force distribution and constant-circulation flight). Variation in
wing span is a fundamental method of generating thrust under the constraint of
constant circulation (Spedding and DeLaurier,
1995⇓); bats, including *P.
poliocephalus*, experience significant variations in tip-to-tip wing span
during steady flight which makes constant-circulation flight plausible (see
also Fig. 3). We also took care
to ensure that the directions of the segmental aerodynamic forces over one
entire wingbeat cycle were such that mean lift equaled weight and mean thrust
equaled total drag. We assume that the small radius of curvature of the
leading edge of the wing produces low leading-edge suction efficiency
(DeLaurier, 1993b⇓); the
consequence of this assumption is that the aerodynamic force must act almost
normal to the wing instead of in a vertical direction.

### Drag

We modeled three distinct drag components that do not arise from
circulation about a wing segment: skin friction over the surface of the wing,
form drag of the body and induced drag from wingtip vortices that exist
because of the finite wing span. The first two drag components are
proportional to the square of the flight speed, whereas the induced drag is
inversely proportional to the square of the flight speed; hence, the induced
drag is an order of magnitude larger than skin friction or form drag on
account of the relatively low flight speed. The maximum segmental velocity can
reach several times the forward flight velocity, and this argues against using
the forward flight velocity in the model; however, we found that the skin
friction drag is negligible for this bat species at the moderate flight speeds
employed in this model and that a more precise treatment is therefore
unwarranted here. As suggested by Norberg (Norberg,
1972a⇓), we assume fully
turbulent boundary layer flow over the entire bat wing because of the small
leading edge radius and the presence of hair near the leading edge of the
wing. The turbulent drag coefficient *C*_{D,f} of a flat, rough
plate is approximately:

(White,
1991⇓), where the mean wing
chord, *c*, is equal to 0.151 m and *l*, the approximate length
of the hair on the bat wing, is 4×10^{-5} m. The drag
coefficient is relatively insensitive to the ratio *c/l*. We write the
total skin friction on both sides of both wings **D**_{f} as:

where *A* is
the mean area of a single wing. Skin friction acts in the plane of each wing
segment and opposes forward flight. We use the mean wing area rather than the
instantaneous wing area because this drag component is too small in magnitude
to warrant an exact treatment. We also assumed that the flight speed
*U* was constant over a wingbeat cycle when evaluating drag; we will
show below that the global acceleration *a*_{sz.q} is too small
to cause significant changes in flight speed.

We approximated body shape as a sphere and assumed that air flow around the
body has a turbulent boundary layer and a typical drag coefficient
*C*_{D,b}≈0.5 (White,
1991⇓). In contrast to the
bodies of birds, those of bats are not streamlined, and can, indeed, resemble
a sphere to some extent because of the proportions of the ribcage. Moreover,
form drag is smaller than the induced drag, which mitigates against searching
for a precise numerical value of the drag coefficient. Our results for this
species (see below) demonstrate that there is not much evolutionary pressure
to drive bat body streamlining from a biomechanical perspective. It follows
that the form drag **D**_{b} is:

where *W* is
the shoulder-to-shoulder width of the body, and π*W*^{2}/4
is the cross-sectional area of the sphere. Form drag acts in the opposite
direction to the direction of flight. Once again, we assumed that the flight
speed *U* is constant over a wingbeat cycle.

Induced drag represents the continuous conversion of freestream momentum
into the wasted angular momentum of wingtip vortices. We assumed that a
constant and uniform induced drag exists across the entire wing span,
including behind the body. However, since we do not consider segmental angles
of attack in our model, we added induced drag explicitly to each wing segment,
rather than implicitly through a downwash angle that alters the orientation of
the aerodynamic force (DeLaurier,
1993b⇓). We calculated the
induced drag on the basis of the mean wing span and an assumed elliptical
distribution of lift along the wing span (Norberg,
1990⇓). The induced drag
**D**_{i} becomes:

where
*M*** g** is total weight,

*k*=1.2 is a constant proposed by Pennycuick (Pennycuick, 1989⇓) and

*S*is the mean wing span. Induced drag acts in the opposite direction to the direction of flight. If the instantaneous wing span is defined as the shortest distance from one wingtip to the other along the wing surface, then one standard deviation in the wing span in 7.5% of the mean.

The total drag **D**_{t} on the bat in forward flight is the sum
of skin friction drag, body form drag and induced drag:
13
where *D*_{t} is the magnitude of the total drag. We have
assumed that the total drag is constant over time because it is an order of
magnitude smaller than the aerodynamic force. On the basis of the data in
Table 1, the ratio of mean
total lift *M*** g** to total drag is

*M*

**/**

*g**D*

_{t}≈9 for this species. The magnitude of the aerodynamic force is found by adding the squares of the lift and drag components, thereby making the lift contribution roughly 80 times more important than the drag contribution.

We assumed for all drag components that the forward flight velocity
*U* is constant over one wingbeat cycle. To evaluate this
approximation, we considered the deceleration caused by the total drag over
one-quarter of a wingbeat cycle. If the deceleration is uniform and
approximately equal to
*a*_{sz}≈*D*_{t}/*M* over a duration
*T*/4, then we can write the fractional change in flight velocity as
approximately:

This value is
approximately -1.4% for *Pteropus poliocephalus* and is 100 times
smaller than *U*, justifying the approximation of constant
*U*.

### Thrust

Bat wings must generate thrust from wing segment circulation to counteract
the drag, as described above. Thrust can be achieved during part of a wingbeat
cycle by temporarily orienting a component of the aerodynamic force in the
direction of flight and always arises from an asymmetry in the wingbeat cycle
(Spedding and DeLaurier,
1995⇓). Sustained horizontal
flight requires that the mean total thrust over one wingbeat cycle equals the
total drag. In the absence of more specific information on thrust generation,
we modeled the variation in net thrust over one wingbeat cycle by making the
following plausible assumptions: the body does not generate thrust, and the
skin friction of a wing segment is overcome by thrust from the same segment.
Therefore, the mean thrust **T**_{f.p} required over one wingbeat
cycle to overcome skin friction on both sides of a wing segment *p*
is:

where the negative
sign indicates a force in the cranial direction and *w*_{p} is
the mean segmental width. Once again, we note that the aerodynamic force is
dominated by lift to the point where detailed variations in drag and thrust
are unimportant in comparison. Consequently, we will define mean wing
parameters that will facilitate a simple and consistent mathematical
approximation of the thrust force. We define *w*_{p}, the mean
width of a given segment *p* over the entire wingbeat, and the mean
wing segmental area *A*_{p} as:

16

17

to be consistent with the
value of skin friction drag above. We calculate that *A*=0.0633
m^{2} from the mean distances between wing landmarks and the geometry
of segmental chords (Table 2).
The difference between the measured wing area
(Table 1) and the value
calculated here is approximately 10%.

We also assumed that form drag and induced drag are each counteracted by
thrust generated uniformly over the mean length of the two wings, 2*B*.
We write the mean thrust needed to overcome form drag and induced drag on a
wing segment *p* as:

where both quantities are
averaged over one wingbeat cycle. We define mean wing length, *B*
as:

20
to recover the values
of the form drag and induced drag given in the previous section. We calculated
*B*=0.431 m from the mean distances between wing landmarks. Since we
distributed these two components of thrust uniformly over the wing span,
whereas lift tapers near the wingtips, the relative importance of net thrust
to lift increases distally.

Passive aeroelastic wings supported by leading edge spars tend to produce
maximum thrust during the middle of the downstroke as a result of wing
twisting that rotates the aerodynamic force towards the direction of flight
(DeLaurier, 1993c⇓). Therefore,
we assumed sinusoidal variation in thrust over the wingbeat cycle, with the
maximum thrust occurring during the middle of the downstroke. We calculated
the magnitude of the thrust explicitly to find the orientation of the
aerodynamic force. In the absence of specific thrust observations, we verified
*a posteriori* that sinusoidal thrust variation does not strongly
affect our simulation results. Ideally, a more precise thrusting function is
desirable, particularly one that captures the 65% downstroke *versus*
35% upstroke cycle, but we did not feel that such an effort was germane to
this work at present. We write the magnitude of the thrust at wing segment
*p* and time *q* as:
21
where f=1/*T* is the flapping frequency, the argument within the cosine
function is expressed in radians, and the time *t*_{q}
associated with the integer *q* is written:
22
where *t*_{mid} is the time at the middle of the downstroke,
and Δ*t*=*T*/40 is the time step. The maximum thrust
2(*T*_{f.p}+*T*_{b.p}+*T*_{i.p})
occurs when Δ*tq*≈*t*_{mid}, and thrust is
equal to zero when Δ
*tq*≈*t*_{mid}+(*T*/2) during the
upstroke. The mean thrust generated by both wings over one wingbeat cycle is
equal to the total drag since the contribution from the cosine function
averages to zero.

### Aerodynamic force distribution and constant-circulation flight

The magnitude of the total aerodynamic force *F*_{a} needed
to sustain forward flight is a combination of weight and total drag components
with a mean value over one wingbeat cycle of:

23

There are three constraints for all plausible spanwise distributions of the total aerodynamic force. First, the increase in circulation proximal to the wingtip is dictated by the mathematical singularity of the wingtip vortex itself: all plausible spanwise distributions of the aerodynamic force are similar near the wingtip, with spanwise circulation becoming zero. Second, flow visualization around a gliding model bat indicates that the body and uropatagium are capable of producing lift (P. Watts, unpublished data). As a result, we postulate the spanwise distribution of aerodynamic force generated by the armwings and the body to be a relatively constant plateau. Third, the constraint that mean lift over one wingbeat cycle equals body weight bounds the possible numerical values of a constant aerodynamic force plateau. An elliptical distribution of the total aerodynamic force over the wing span is a common and convenient choice for engineers designing subsonic fixed-wing aircraft and even ornithopters (DeLaurier, 1993b⇓). It produces minimum induced drag as well as a uniform downwash across the wing span; thus, an elliptical distribution of aerodynamic force over the wing span would be energetically favorable for flying mammals. Given the three constraints mentioned above, an elliptical distribution of aerodynamic force is an appropriate and plausible way of connecting an aerodynamic force plateau to wingtip vortices. To balance weight, any other plausible distribution that is less than an elliptical distribution somewhere along the span must also be greater than the elliptical distribution elsewhere along the span. We expect errors incurred by actual deviations in the spanwise distribution of aerodynamic force from the chosen elliptical form to be small in magnitude compared with other approximations made in the model.

Spedding (Spedding, 1987⇓)
used stereoscopy of small helium bubbles in the wake of a kestrel, *Falco
tinnunculus*, to reveal that flapping flight need not be accompanied by
the shedding of spanwise vortices from the trailing edge of a wing. The
absence of vortex shedding implies a constant circulation about each spanwise
segment of the wing in fast forward flight. Consequently, constant circulation
implies that each spanwise wing segment maintains an aerodynamic force that
changes orientation but not absolute magnitude during the wingbeat cycle. It
has been demonstrated that bats can use a constant-circulation gait during
fast forward flight and that bat wings possess sufficient passive and active
control to achieve constant-circulation flight over an entire wingbeat cycle
in fast forward flight (Rayner et al.,
1986⇓). The approximately 6 m
s^{-1} flights from which our data were collected may not fall
strictly within the fast flight category; Norberg and Rayner (Norberg and
Rayner, 1987⇓) estimate that
the maximum range and minimum power flight speeds for flying foxes of this
mass are approximately 10 and 7.5 m s^{-1}, respectively, on the basis
of the general relationship they estimated between body mass and flight
velocity. In any case, neither the constant-circulation theory of flapping
flight nor the flow visualization on which it is based is capable of providing
the spanwise distribution of the aerodynamic force. Therefore, an elliptical
aerodynamic force distribution combined with the constant-circulation
hypothesis appears to be a reasonable starting point for modeling
purposes.

In our model, we defined the magnitude of the aerodynamic force on a wing
segment *p* as:
24
where Γ_{pq} is the circulation about the midline of the wing
segment and *w*_{pq} is the instantaneous width of the wing
segment. We invoked the constant-circulation hypothesis when we made the
magnitude of the aerodynamic force Γ_{pq}*w*_{pq}
constant over time for each wing segment. There are two possible
interpretations of constant circulation. The global constant-circulation
hypothesis requires that the sum of circulation along the entire wing span
remain constant during flapping flight. The local constant-circulation
hypothesis requires that the sum of circulation between two close material
markers along the wing span remains constant. The sum of circulation along the
wing span is directly proportional to the aerodynamic force. The global
hypothesis allows for a spanwise redistribution of circulation, which would
lead to visible vortex rings being shed from the trailing edge. An absence of
such vortex rings implies the stronger local hypothesis, which we are using
here. Because we are using sums over wing segments instead of integrals along
wing span in this work, we make our material markers coincident with the wing
segment boundaries and enforce a constant aerodynamic force over each wing
segment. Because a bat wing undergoes significant folding in flight, we have
chosen the middle of the downstroke as the most appropriate wing configuration
at which to distribute the aerodynamic force along the wing span. Of the 40
time steps simulated by the model, mid-downstroke is identified as the integer
*q*=*Q*, where Δ
*t*_{q}≈*t*_{mid}=00.38*T* or
*Q*=15.

We approximated the smooth elliptical distribution traversing each wing
segment with the value of the circulation at the segment midline, resulting in
a spanwise distribution of circulation that resembles a descending staircase
over the entire wing length from root to tip
(Fig. 4). The circulation Γ
_{pq} is calculated according to the *x* axis position of
the midline of each wing segment relative to the center of the body at the
instant of the middle of the downstroke. Because the maximum circulation
exists midway between the shoulders, the circulation about a wing segment
*p* is given by:

(Norberg,
1990⇓), where
*x*_{pQ} is the position of the midline of a wing segment
relative to the shoulder, *B*_{mid} is the instantaneous wing
length at *t*_{mid}, *F*_{a} is the total
aerodynamic force and *C* is a multiplicative constant. This equation
is strictly valid only for straight wings with constant wing span. Since we
have made several approximations while distributing the aerodynamic force
along the wing span at one instant in time, and since the wing moves and
changes shape during a wingbeat, circulation must be multiplied by an
aerodynamic constant *C* to ensure that the mean lift exactly balances
body weight over one wingbeat cycle. We solved for the aerodynamic constant
*C* by iteration. We approximated the circulation about the body (from
shoulder to shoulder) by:

the maximum
circulation in the elliptical distribution at *x*=-*W*/2.
Substituting this value of the circulation into the definition of the
aerodynamic force generated by the body yields:

where *W* is
body width and *F*_{a} is the total aerodynamic force, a simple
function of total drag and body weight.

### Net drag and lift

We defined net drag by adding the drag and thrust acting on each wing
segment. The net drag on wing segment *p* at time *q* is given
by:
28
and acts in the plane of the wing segment opposing the direction of flight. We
calculated the magnitude of the lift from the mathematical relationship
between the magnitude of the aerodynamic force and the magnitude of its two
orthogonal components, lift and thrust. Hence, the lift on wing segment
*p* at time *q* is written:

where the lift
component *F*_{l.pq} of the aerodynamic force acts normal to
the plane of a wing segment. We assumed that the lift on each segment acts at
a fixed proportion, *P*_{c}≈4, of the chord,
*c*_{p}, behind the leading edge of a wing segment. This is a
typical location for the center of lift of a thin airfoil with a significant
parabolic camber and angles of attack between 5 and 15° (Katz and Plotkin,
1991⇓).

### Internal force

The internal force is the force carried by the wing structures that enables
the wing to undergo the observed accelerations and to resist the external
forces applied to the wing. We invoke the common form of Newton's second law
of motion for wing segment *p* at time *q* so that:
30
where the internal force **F**_{pq} is not given any subscript. We
did not calculate the internal force explicitly in the model since we were
interested in the relative contribution of inertial and external forces to
bone stresses. Moreover, the inertial and external forces act at different
locations on a given wing segment, requiring both a net internal force and a
net internal moment about the center of mass for each wing segment. Instead,
we relied on the definition of the internal force as the inertial force minus
the gravitational force, added mass force, lift and net drag, and computed the
contributions from these forces separately. We also calculated the component
of the internal force associated with tension in the plagiopatagium (see
below, section on plagiopatagium tension).

### Global acceleration

To calculate the body's global acceleration, we employed the estimate of total internal force transferred from the wings to the body. We assumed that the position of the shoulder joint relative to the center of mass is fixed then, working from the inertial reference frame of the global axes, approximated the force transfer from the wings to the body by summing the segmental internal forces from the wing root to the wingtip at each instant of time. The components of the total internal force transferred from both wings are:

where the *x* axis
components from both wings are equal and opposite, resulting in no
acceleration of the center of mass along the *x* axis at any time in
forward flight. To complete the force balance, we calculated the mass of the
body minus the wings *m*_{b} as:

The vertical force
balance on the body supported by both wings at time *q* is:
34
where *a*_{sy.q} is the vertical component of the global
acceleration, and *F*_{a.0} is the aerodynamic force generated
by the body. The horizontal force balance provides at time *q*:
35
where *a*_{sz.q} is the horizontal component of the global
acceleration. Form and induced drag of the body have already been accounted
for in the net drag of each wing segment. We solved these two equations for
the two components of the global acceleration iteratively because the
equations for *F*_{y.pq} and *F*_{z.pq} contain
the global acceleration components within inertial terms.

### Level flight criterion

As mentioned above, we ensured that the mean lift was equal to body weight
over the wingbeat cycle by adjusting the aerodynamic constant *C*.
Without explicit knowledge of the correct value of *C*, we adopted the
following approach: we made two initial estimates near unity, observed the
resulting trend in the mean vertical component of the global acceleration and
used these results to guide more accurate estimates of the aerodynamic
constant *C*. To do this, we averaged the equation for the vertical
component of the global acceleration over one wingbeat cycle to obtain:

where the correct
value of the aerodynamic constant *C* makes the left side of this
equation equal to zero. Under this condition, lift from the wings and body
balances body weight over one wingbeat cycle. We used a Newton—Raphson
iterative convergence scheme to compute values for *C* until the left
side was less than 10^{-10} at each iteration. The
Newton—Raphson scheme converged in less than 60 iterations to final
values of *C*≈1.2 as well as final values of global accelerations.
The aerodynamic constant is greater than unity largely to compensate for the
significant wing folding that occurs at the end of the downstroke.

### Plagiopatagium tension

Lift acting on the armwing causes the plagiopatagium to billow and thereby increases the tension of the skin. Because of this coupling between lift and skin tension, we replaced total armwing forces with skin internal forces; skin force represents the sum of skin tension integrated along the chordwise direction and does not include the frictional drag, which is accounted for separately through the net drag. We therefore developed a simple model of passive plagiopatagium deflection due to quasi-steady lift to determine the internal forces of the skin that are transferred directly to, and indirectly through, the humerus and radius.

To calculate the magnitude and location of this force, we assumed that the
skin of the plagiopatagium is linearly elastic and orthotropic, that skin
thickness remains constant during the wingbeat cycle, that shear stresses in
the skin may be neglected and that the intrinsic musculature of the membrane,
particularly the mm. plagiopatagiales, does not affect skin tension near the
fifth digit during the wingbeat. This final assumption comes from our
hypothesis that the muscles will only have a localized effect on membrane
tension, perhaps by influencing local camber for aerodynamic purposes or by
damping oscillations in the wing membrane *via* modulating skin
stiffness. These localized effects will probably have little influence on skin
tension measured elsewhere in the wing. The nature of elastic problems is such
that a local perturbation of the global solution is rapidly smoothed out of
existence. This local *versus* global dichotomy is a consequence of the
elastic partial differential equation being elliptical: the stress at every
point on the wing membrane is in some way an average of the stresses
surrounding that point. Therefore, elastic materials quickly erase any
evidence of local perturbations.

We characterized the plagiopatagium as a rectangular sheet coplanar with
the plane connecting the shoulder, elbow and carpus if the wing membrane was
not deformed by lift (Fig. 5).
The long axis of this rectangular membrane is parallel to the line connecting
the shoulder to the carpus. The rectangle's chordwise dimension is the mean
armwing chord *c*_{a} and its spanwise dimension is the
instantaneous shoulder-to-carpus distance *b*_{a}. For this
analysis, we used a coordinate system with its origin at the carpus, its
*x* axis parallel with the long edge of the rectangle and positive
proximally, and its *z* axis parallel with the fifth digit and positive
caudally (Fig. 5). The location
and orientation of the total skin internal force *F*_{s}
justifies the approximation of the plagiopatagium as a rectangle since the
line of force action passes distal of the armwing bones (see below). Newton's
third law requires that the plagiopatagium induce a force normal to the
armwing that is equal to the lift that it is replacing. Part of the lift is
accommodated by the angle tan^{-1}*s* that the force
*F*_{s} makes with the plane of the fifth digit and the body,
and we distributed the remainder of the lift evenly along the armwing bones
acting normal to the armwing.

We estimated the spanwise (*E*_{xx}) and chordwise
(*E*_{zz}) elastic moduli of *P. poliocephalus*
plagiopatagia from published mean values of moduli of plagiopatagia of other
bats, yielding *E*_{xx} =3 MPa and *E*_{xzz}=37
MPa; these values correspond closely to the elastic moduli found for the
plagiopatagia of *A. jamaicensis*, the only large fruit-eating bat from
which these data have been collected (Swartz et al.,
1996⇓). However, because this
species is an order of magnitude smaller in body mass than *P.
poliocephalus*, we also explored the effects of lower moduli on force
transfer. We further assumed the Poisson's ratio ν_{zx}=1 as a
typical value for skin composed of crossed fiber sheaths (Frolich et al.,
1994⇓). Symmetry of the skin
stress tensor associated with mechanical equilibrium requires that ν
_{xz}=(*E*_{xx}ν_{zx})/*E*_{zz}=0.8,
where ν_{xz} is Poisson's ratio of the plagiopatagium due to
loading along the *x* axis.

Although inertial forces exerted on the skin can probably approach the
magnitude of aerodynamic forces near the carpus during periodic rapid wing
acceleration, we modeled only the skin deformation induced by the
instantaneous lift acting on the plagiopatagium. We divided total armwing lift
by the instantaneous plagiopatagium area
*b*_{a}*c*_{a} to find mean pressure differences Δ
*P* (proportional to the wing loading), assumed to act uniformly
at each instant of time over the plagiopatagium as the source of membrane
deflection. We adopted a quasisteady model of membrane deflection out of the
plane by neglecting the retarding effect of membrane inertia on changes in
plagiopatagium deflection over time. We computed the linear solution of the
partial differential equation governing membrane deflection *Y* using
the method of separation of variables by assuming uniform values of the skin
tensions *T*_{x} and *T*_{z} in the spanwise and
chordwise directions respectively. To capture the most important contribution
to deflection (and hence skin tension) and yet maintain a tractable analytical
solution, we approximated the complete linear solution of the partial
differential equation by the first eigenmode solution corresponding to the
boundary conditions that there is no membrane deflection along the distal,
cranial and proximal edges, and that the slope of the membrane deflection be
zero along the caudal edge of the rectangle. The linearized partial
differential equation applies whenever lateral membrane deflections are small
enough to be considered remaining in the undeflected plane, which is an
assumption that must be verified *a posteriori*. We find the
approximate solution for membrane deflection out of the plane:

where the skin
tensions *T*_{x} and *T*_{z} are functions of
*Y* through the orthotropic constitutive relationships:
38
and
39
where we assumed uniform plane stress across the membrane thickness
*t*_{s}=0.2 mm. We evaluated the spanwise (ϵ_{xx})
and chordwise (ϵ_{zz}) skin strains by assuming that, when under
no tension, the membrane had a chordwise length *c*_{a} and a
spanwise length *b*_{min}, the minimum value of
*b*_{a}. It follows from the definition of an arc length that
the strain ϵ_{xx} is:

and that membrane
deflection out of the plane both determines and is determined by skin tension.
We do not provide a sign convention for Δ*P* or *Y*
because the square of the derivative ∂*Y*/∂*x* makes
these quantities positive. We did not calculate the strain ϵ_{zz}
since we were interested in the tension *T*_{x} along the edge
*x*=0, where the strain ϵ_{zz} is identical to zero and
*T*_{z}=*T*_{x} from the definitions of the
Poisson's ratios.

We calculated the skin tension *T*_{x} at 20 discrete
locations along the distal edge of the plagiopatagium using a
Newton—Raphson convergence scheme at each location to solve iteratively
the balance between deflection and tension. We found uniform skin tensions
along the fifth digit except for a brief period during the upstroke at which
time the skin tension approached zero and also became non-uniform. To
calculate a particular *T*_{x}, we transformed the integral for ϵ
_{zz} into a complete elliptical integral of the second kind,
*E(m,π*/2):

where
*m*≡*s*^{2}/(1+*s*) and *s* is the
value of the derivative ∂*Y*/∂*x* evaluated at
*x*=0. We estimated *E*(*m*,π/2) from a fifth-order
polynomial curve fit that had a standard deviation of 0.12% about the exact
integral values over 0≤*m*≤1 (Abramowitz and Stegun,
1965⇓). We found that m≈0.01
and *E*(0.01,π/2)≈1.57 for most of the wingbeat, and therefore
employed a linear solution since tan^{-1}≈5° is the angle that
the membrane makes with the undeflected plane along the edge *x*=0. The
total skin internal force *F*_{s} pulling proximally on the
fifth digit is the integral of the tension *T*_{x} along the
*z* axis from *z*=0 to *x*=*c*_{a} and can
be approximated by
*F*_{s}=*T*_{x}*c*_{a} over most
of a wingbeat (Fig. 5). Since
the fifth digit is free to rotate about the carpus, we also calculated the
total moment *M*_{s} induced by the skin force
*F*_{s}. The correct point of action of the total skin internal
force is found from
*z*=(*M*_{s}/*F*_{s})≈0.5*c*_{a}
in this special coordinate system.

### Forces, moments and stresses on the wing skeleton

To analyze the effect of flight-related forces on the skeleton, we
designated the point at which an inertial or external or skin membrane force
vector **F**_{f} acts on a given wing segment as
(*x*_{f}, *y*_{f}, *z*_{f}) and
the point on the skeleton about which we seek the moment **M**_{0}
as (*x*_{0}, *y*_{0}, *z*_{0}).
By varying the location of the point for which the analysis is carried out, we
may then analyze the skeletal loading and bending at any location on the wing
skeleton of the armwing that is biologically significant and through which a
known fraction (usually taken as 100%) of the remaining internal force is
transmitted. Although the analysis applies generally to any location along the
armwing, we focused here on the moments about the humeral and the radial
midshafts, to carry out comparisons of model estimates with empirically
measured values (Swartz et al.,
1992⇓). We defined three moment
arms:
42
43
44
that relate the distance between the applied force and the location of the
moment. We calculated the moment **M**_{0} about the point
(*x*_{0}, *y*_{0}, *z*_{0})
induced by the force **F**_{f} at (*x*_{f},
*y*_{f}, *z*_{f}) as:
45
46
47
(Meriam and Kraige, 1987⇓). To
calculate the total force and moment vectors at bone midshafts, we summed the
inertial and external forces **F**_{pq} and their respective
moments *M*_{0.pq} over all wing segments and fractions thereof
that are distal to the bone midshaft using the appropriate locations
(*x*_{f}, *y*_{f}, *z*_{f}) for
all five forces within a given wing segment. We then subtracted the internal
force and moment about the bone midshaft associated with the skin force since
it is transmitted directly to the body. We note that our calculations
implicitly contain close approximations of the moments of inertia of the wing
through the summation of segmental inertial forces, although we neglected the
small moments of inertia induced by the finite width of the wing segments
(Thollesson and Norberg,
1991⇓).

From the total force and moment vectors at a bone midshaft, we calculated
stresses at dorsal and ventral sites on a bone's circumference treating the
bone as a beam of circular cross section (Popov,
1978⇓). We considered the
center of the bone cross section as the origin of a local coordinate system
parallel to the global coordinate system
(Fig. 2,
Fig. 6). We converted total
force and moment components written in the local *x, y, z* coordinate
system into a primed *x′, y′, z′* coordinate system
oriented to the bone long axes. This primed system shares the origin of the
local coordinate system, but possess an *x′* axis oriented
parallel to the bones' long axes (positive distally), a *z′* axis
in the plane of the wing (positive caudally) and a *y′* axis
normal to the wing plane (positive dorsally)
(Fig. 6). We defined the angle λ
as the angle between the armwing and the horizontal *x* axis.
The angle λ defines a new *y′* axis normal to the armwing.
We also defined ψ as the angle by which we rotate the *y′*
axis to align the *x′* axis with the long axis of the humerus or
radius. The coordinate transformations from the local axes to the primed axes
are:
48
49
50
The same transformation is used to relate the primed components
*F*_{x′}, *F*_{y′},
*F*_{z′} and *M*_{x′},
*M*_{y′}, *M*_{z′} to their
respective unprimed components; we reapplied the transformation at each
instant of time since the angles λ and ψ change throughout the
wingbeat cycle. To calculate the angle ψ, we applied the law of cosines to
the triangle formed by humeral length *L*_{h}, radial length
*L*_{r} and the instantaneous distance from the shoulder to the
carpus. This yielded the angles interior to this triangle, which then provided
the position of the elbow joint and the angle ψ for either the humerus or
the radius.

We assumed that the humerus and radius are straight circular cylinders of
inner radius *R*_{i}, outer radius *R*_{o} and
length *L*, and that the medullary cavity has no capacity to resist
bending or torsion. We also assumed that the neutral axis intersects the
centroid of the cylinder's cross section. To calculate stresses, we balanced
the total force and moment vectors at midshaft with longitudinal normal
stresses σ_{L} and shear stresses τ_{1,2}
distributed within the bone at the surface of a perpendicular midshaft cut.
When integrated over the proximal surface of the cut, the normal and shear
stress distributions represent three force and three moment components that
are equal in magnitude but opposite in sign to the total force and moment
components. We considered tensile normal stresses and left-handed (or
clockwise) shear stresses as positive (Fig.
7). We also assumed that the resulting stresses can be
superimposed linearly given that superposition of stresses is quite accurate
up to approximately 1% strain; we note that the peak humeral and radial
strains in *P. poliocephalus* rarely exceed 0.3% (3000 μϵ)
(Swartz et al., 1992⇓).

Assuming that a bending moment induces a linear increase in normal stress σ
_{L} with increasing distance from the neutral axis, we
calculated the magnitude of the normal stress induced by
*M*_{y′} as:

where *J*, the
polar second moment of area of the bone cross section, is π
(*R*_{o}^{4}-*R*_{i}^{4}),
and we are free to choose
*R*_{i}≤*z′*≤*R*_{o}. For
positive *M*_{y′}, the cranial face of the bone is in
tension while the caudal face is in compression, as indicated by the sign of σ
_{L}. The magnitude of the normal stress induced by
*M*_{z′} is:

where positive
*M*_{z′} indicates tension along the bone's dorsal
surface. Assuming that stress increases linearly from the *x*′
axis in the radial direction, the magnitude of the shear stress induced by the
torsional moment *M*_{x′} is:

where
*R*_{i}≤*R*≤*R*_{o}. This shear
stress is uniform around the bone circumference, and the maximum values of
both the normal and shear stresses occur at the bone outer radius
*R*=*R*_{o}.

The normal stress induced by the force *F*_{x} was
calculated by assuming that the force is distributed evenly over the cross
section, giving rise to a normal stress:

where
*A*_{c} is the bone's cross-sectional area. Forces acting
perpendicular to the long axis of the bone give rise to shear stresses that
are not distributed evenly along the bone cross section (Beer and Johnston,
1981⇓). The magnitude of the
shear stress at *y*′=0, *z*′=±R_{o}
due to the internal force *F*_{y′} is:

where, for positive
*F*_{y′}, shear stress is positive on the cranial
(leading) and negative on the caudal (trailing) surface of the bone. The
magnitude of the shear stress due to the internal force
*F*_{z′} at *z*′=0,
*y*′=±*R*_{o} is:

where, for positive
*F*_{z′}, shear stress is positive on the dorsal and
negative on the ventral surface of the bone. The two magnitudes of τ
_{1,2} given here apply to only two points, each along the outer
perimeter of the bone.

Swartz et al. (Swartz et al.,
1992⇓) assumed that compact
cortical bone is a transversely orthotropic material and used published values
for compact bone moduli and Poisson's ratio (Carter,
1978⇓) to convert empirical
measured strains into estimated stresses. They used a value of
*E*_{11}≈15 Gpa, where *E*_{11} is the
compressional longitudinal elastic modulus of bone, a typical value for
mammalian compact cortical bone (Currey,
1984⇓). Since there are no
obvious sources of transverse stress σ_{T} applied to the
humerus or radius, we assumed that applied σ_{T}=0. Hence, the
longitudinal strain is simply ϵ
_{1}=σ_{1}/*E*_{11} and transverse
strains ϵ_{2}=-0.46ϵ_{1} exist due to Poisson's
ratio effects alone, even in the absence of externally applied transverse
forces.

## Results

### Model predictions

#### Humeral and radial stresses

We used the model described here to estimate stresses at the mid-dorsal and mid-ventral locations on both the humerus and the radius, and at the cranial (leading edge) and caudal (trailing edge) midshaft of the humerus; we did not calculate cranial and caudal stresses on the radius because we did not have empirical data with which we could compare the model results.

During the portions of the wingbeat cycle in which the wing is extended, the wing membrane is placed under significant tension. This tension exerts forces on the humerus and radius that tend to bend them in the plane of the wing in a cranially convex manner, placing the cranial bone surfaces in tension and the caudal surfaces in compression. Values of model predictions for stresses on the cranial and caudal surfaces of the bones are therefore strongly affected by input values of skin moduli (Fig. 8A; see also sensitivity analysis, below). Because this bending is largely restricted to the plane of the membrane, skin modulus input has little effect on estimates of shear and longitudinal dorsal and ventral normal stresses (Fig. 8B,C).

To date, no studies of the mechanical properties of bat wing membranes have
included large megachiropterans (Studier,
1972⇓; Swartz et al.,
1996⇓). Therefore, we began our
analyses using values for mean spanwise and chordwise skin moduli measured
from bats an order of magnitude smaller in body mass than our model species.
On the basis of a physical inspection of the wing membranes of many bat
species, we expected these values to be considerably higher than the true
moduli of *P. poliocephalus* wing membrane skin, and carried out
analyses using lower values for moduli as well, maintaining the empirically
measured ratio of spanwise to chordwise modulus
(Fig. 8). When the model
employs a spanwise modulus of 0.5 MPa, one-sixth of the mean spanwise
stiffness of the wing membranes of smaller bats (Swartz et al.,
1996⇓), the magnitudes of
estimated bone stresses on the cranial and caudal bone surfaces are most
similar to those on the dorsal and ventral surfaces, so we conducted all
further analyses with this lower value for the skin moduli; this is reflected
in the results for the cranial and caudal surfaces of the humerus.

Predicted stresses at all bone sites except the cranial humerus share an overall pattern of change with respect to the wingbeat cycle: stresses reach a peak just after the initiation of downstroke and another peak of similar, slightly greater or slightly lower magnitude at the middle of downstroke; stresses then decrease to their lowest values shortly after the initiation of upstroke (Fig. 9). On the cranial aspect of the humerus, the longitudinal stress peaks just after the beginning of downstroke, exhibits a plateau throughout most of the remainder of the downstroke and then decreases steeply at the initiation of upstroke in concert with the stress minima of the other sites. Longitudinal stresses at all sites are higher than shear stresses throughout the greater part of the wingbeat cycle with the exception of the downstroke—upstroke transition, at which time longitudinal stress values approach zero.

The underlying basis of the similarity in the stress profiles at all anatomical locations modeled becomes clear when total stress is partitioned into components due to gravitational, inertial, added mass, net drag, handwing lift and armwing membrane (plagiopatagium) tension stresses (Fig. 10, Fig. 11). The inertial stresses display the pattern of timing observed for net stresses at all sites other than the cranial humerus, displaying peaks just after the initiation of downstroke and at mid-downstroke, and minima at the downstroke—upstroke transition. Although the general shape of the inertial stress curve for the cranial humerus is similar in overall form to that of the dorsal and ventral bone surfaces, the magnitude of inertial stress (+5 to -10 MPa) is far smaller here than stress due to tension in the plagiopatagium, so plagiopatagium stresses dominate. Skeletal stresses due to gravity, added mass and drag are virtually always very low, typically below 4 MPa and often close to 1 MPa. The only exception to this pattern is in the longitudinal stress at the cranial mid-humerus; in this case, the thrust engendered at mid-downstroke is large enough to generate significant bending in the plane of the wing, decreasing tension on the cranial surface by up to approximately 7 MPa or 20 % of total stress at this location.

Stress due to plagiopatagium tension is also relatively small, generally around 2-5 MPa but reaching contributions as high as 10 MPa to the total longitudinal stress on the dorsal and ventral surfaces of the humerus. The contribution of plagiopatagium stress to the total stress at the cranial edge of the humerus is much greater, as high as 40 MPa, for the geometric reason noted above. If the model used a higher value for spanwise membrane modulus, this contribution would be even higher. Stresses in the humerus and radius due to the lift generated by the handwing are significant (5-20 MPa) and show a consistent plateau with little variation in magnitude throughout the downstroke.

### Global body motion

The model predicts that the bat's center of mass has a positive (upward)
vertical acceleration through the greater part of the downstroke, downward
acceleration through the end of the downstroke and beginning of the upstroke
and upward acceleration once more during the second half of the upstroke
(Fig. 12). Downward
acceleration reaches a maximum of nearly 13 N; given that this animal's mass
is approximately 800 g, this acceleration is close to 2 ** g**.
This large acceleration occurs because the wing is folded such that some
aerodynamic force acts downwards near the top of the upstroke. In addition,
the upward wing acceleration near the end of the downstroke causes further
downward acceleration of the body.

During the middle half of the downstroke, the body is accelerating forward, but it loses its forward acceleration and gradually decelerates during the end of the downstroke and first half of the upstroke (Fig. 12). Near the upstroke—downstroke transition, there is a brief forward acceleration followed by a brief deceleration before the acceleratory phase of the downstroke begins. The magnitude of the fore—aft component of the acceleration of the center of mass (±2 N) is small compared with the vertical acceleration. As a consequence, forward flight speed varies little during single wingbeats in level flapping flight.

### Model validation: comparisons with empirical data

#### Humeral and radial stresses

To determine the degree to which our model accurately estimates bone
stresses, we compared midshaft longitudinal and shear stresses calculated from
the model with values computed from directly measured *in vivo* wing
bone strains of bats from the same study population. In this earlier study
(Swartz et al., 1992⇓), maximum
and minimum principal strain magnitudes and orientations throughout the
wingbeat cycle were computed from data obtained from rosette strain gauges
surgically implanted around the midshaft circumference of the humerus and
radius; post-operative recordings, synchronized with video recordings, were
made as bats flew the length of a 30 m flight cage at moderate speed.

We found that the model can predict effectively not only the general shape of the stress curve, but also the timing of both major and secondary stress peaks in relation to the wingbeat cycle (Fig. 13). Moreover, stress magnitude estimates correspond quite closely to empirically measured values, with predicted values deviating from recorded values by approximately 5-30 % for longitudinal and shear stress peaks.

Deviations between the predicted and recorded values are sometimes greater during the lower-stress portions of the loading cycle and occasionally differed by as much as a factor of 2-3. These deviations arise partly from an `offset' that appears in some of our results: predicted and measured stress plots are very closely matched in shape but offset by 20-25 MPa. This phenomenon is particularly striking at the cranial surface of the humerus (Fig. 13).

Empirically recorded bone strains indicated relatively large transverse
stresses at the humeral and radial midshafts (Swartz et al.,
1992⇓). The constitutive model
for bone presented by Carter (Carter,
1978⇓) is almost
volume-conserving and predicts ϵ_{2}=0.46ϵ_{1} from
the Poisson's ratio alone in the absence of externally applied transverse
stress. It is therefore possible that transverse stresses observed in bat
bones are largely due to Poisson's ratio effects. We computed transverse
stresses from the model as σ_{T}=0.46σ_{L}, where σ
_{T} is transverse and σ_{L} is longitudinal
midshaft bone stress, and then compared these with transverse stress values
calculated from the rosette strain gauge data
(Fig. 14). Correspondence
between the predictions and measured values is excellent and detailed in some
cases, but only moderate in others.

### Global body motion

Because the model predicts the net vertical force exerted on the bat's center of mass during the wingbeat, a comparison of the vertical changes in position of the bat's center of mass calculated from the model with those measured directly from film provides a strong test of model accuracy. Body oscillations are slightly greater in the wind-tunnel flight than in the simulation but, overall, oscillations computed from the simulation and measured directly from film show a high degree of correspondence (Fig. 15).

### Sensitivity analysis

The model's output values for bone stresses and global body motion descriptors are sensitive to various inputs and model assumptions to greatly varying degrees. To assess the importance of these effects, we carried out a sensitivity analysis in which we varied each of 10 model inputs by ±5 % and computed the resulting changes in 24 model outputs (Table 4).

As noted above, the elastic modulus of the wing membrane skin in the spanwise direction has a significant effect on estimates of skeletal stress; a 10 % change in modulus can produce a nearly equal change in bone stress (Table 4; see minimum longitudinal stress, cranial and caudal radius). The sensitivity analysis indicates that this effect arises from the influence of skin modulus on skin tension. Bone shear stresses are also rather sensitive to aerodynamic parameters, e.g. the timing of peak thrust within the downstroke, the timing of mid-downstroke within the entire downstroke, the location of the center of aerodynamic pressure and the added mass coefficient. The humerus and radius are not, however, affected equally; humeral shear stress is more than an order of magnitude more sensitive than radial shear stress. Longitudinal stresses are influenced by the same parameters, but generally to a lower degree. The cranial and caudal maximum longitudinal stresses, however, are moderately sensitive to the location of the center of pressure (Table 4).

In contrast, the model results are quite robust to variation in many model input parameters. Whole-body outputs (accelerations, velocities and displacements of the animal's center of mass) are virtually unaffected by changes in model inputs (Table 4). Longitudinal bone stresses are not influenced by changes in values of skin Poisson's ratio, skin friction or form drag and are only minimally affected by variation in the added mass coefficient.

## Discussion

### Model validation

The flapping flight of vertebrates and the aerodynamics and mechanics
underlying this mode of locomotion are extremely complex phenomena. The
supporting structures of the vertebrate wing are made from composite materials
whose internal architecture far outstrips most engineered materials in
complexity of design. Many of these materials, particularly muscles, tendons,
skin and, in the case of birds, feathers, are not linearly elastic and possess
mechanical properties that can change dynamically through the wingbeat cycle
*via* modulation of tensions in any of a large number of wing muscles.
Wing surfaces move through space in a manner that is rarely spatially simple
and that varies a great deal along an axis from shoulder to wingtip. The
forces experienced by the wing also vary among anatomical regions and
throughout the wingbeat cycle and comprise more numerous and less easily
estimated components than limb forces in typical modes of terrestrial
locomotion. Together, these and other considerations suggest that computer
modeling may be an especially productive approach to gaining greater insight
into how vertebrates in general, and bats in particular, fly.

The greatest strength of computer modeling approaches — their ability to simplify intractably complicated problems and to allow investigators to vary parametrically single constituent elements that cannot be teased apart experimentally — is counterbalanced, however, by their greatest weakness. If the phenomena we model are those whose complexity defies ready empirical study, how are we to assess the validity of our models? Here, we present a model that facilitates analysis of some of the tremendous complexity of bat flight dynamics while retaining a strong link to empirical studies. This model is sufficiently detailed to predict both the net motion of a flying animal's center of mass and the stresses developed in the bones of the wing throughout the wingbeat cycle. As a consequence, we are able to assess the validity of our model with a degree of rigor unusual for computer models of animal flight.

We find that the model described here is able to make excellent estimates of aspects of wing and whole-body mechanics from kinematics and structurally important aspects of wing form and mass distribution. Indeed, a comparison of model results with empirically measured strains and of model-predicted vertical oscillations of the center of mass with those measured indicates that the model is surprisingly realistic (Fig. 13, Fig. 15). The model is able to predict the overall pattern of bone stress change during the wingbeat, the proper modes of stress (tension, compression, shear), some of the small-scale changes in stress during a wingbeat and, with varying precision, the magnitude of the wing bone stresses. Hence, it is reasonable to conclude that the model captures many of the features of horizontal constant-velocity flight that are important to flight performance. The unverified assumptions built into the model are therefore either likely to be confirmed by future direct measurements or are of relatively little importance to the mechanical behavior of the wing in flight.

In addition to longitudinal stresses, the model is also able to estimate
the stresses oriented perpendicular to the long axes of the humerus and radius
from a reasonable estimate of the Poisson's ratio of compact cortical bone
(Fig. 14). The relatively
large transverse stresses observed in the armwing bones of *P.
poliocephalus* are, as a consequence, best understood as arising from this
aspect of bone's mechanical properties rather than from an off-axis pull by
the wing musculature or tension within the wing membrane, as previously
suggested (Swartz et al.,
1992⇓). Indeed, decomposition
of total force exerted on the wing into its constituents demonstrated that the
longitudinal stresses in the humerus and radius due to plagiopatagium tension
are quite low, typically less than 5 MPa, and transverse forces should have
roughly similar values. At some points in the wingbeat cycle for some bones,
however, model-based predictions of transverse stresses are considerably lower
than recorded values. These discrepancies may well indicate instances where
local muscle pull strongly influences recorded strains, and they may be
employed in the future to develop specific hypotheses concerning the activity
of particular wing muscles.

Similarly, the model predictions of vertical movements of the bat's center of mass are closely, although not perfectly, matched by empirical measurements (Fig. 15). The vertical oscillations measured from the film records also display a significant amount of variation. This variation and the lack of precise matching of the two curves is probably due in part to the derivation of the empirical values from the mean of a number of wingbeats whose amplitude and frequency were not identical. Although, to date, there have been no studies of variation in kinematics and/or body motion in bats in relation to wing kinematics, flight velocity, etc. at the level of detail needed to assess such features as wingbeat amplitudes on this scale, we believe that in the future it may be possible to relate kinematics to variation in aspects of flight performance such as the lift generation reflected in the motion of the center of mass.

### Sensitivity analysis

A notable feature of the sensitivity analyses presented here is the dependence of stresses on the cranial and caudal aspects of the bones and the skin tension on the values for skin moduli entered into the model (Fig. 8). Although relatively little is known to date about the mechanical properties of wing membrane skin (Studier, 1972⇓; Swartz et al., 1996⇓), our analyses indicate that they may play an important role in flight mechanics and aerodynamics. Although some previous work has characterized wing membrane mechanical properties across taxa and among wing regions (Swartz et al., 1996⇓), these comparisons were based solely on the linear region of a complex, J-shaped stress—strain curve. The analyses presented here suggest that wing membranes may not reach these relatively high stiffness values during level flight at moderate speed and may operate primarily in the toe region of the curve. Further work will be needed to characterize more accurately this non-linear region of the wing membrane stress—strain curve and to determine the stiffness of the wing membrane during natural flight behaviors.

Beyond this effect of skin modulus, the overall model results are quite robust to variation in a number of input parameters (timing of peak thrust and mid-downstroke, etc.) and extremely robust to a number of inputs (skin Poisson's ratio, coefficient of friction of the skin, form drag) (Table 4). Within this overall pattern, it is striking that the humerus is approximately an order of magnitude more sensitive than the radius to the timing of the mid-downstroke within the wingbeat and the shear stresses due to the location of the center of mass. It is likely that this is a function of a combination of the cross-sectional geometry of the bone and the position of the bone within the wing which, in turn, influence the typical mode of loading of each of the bones. The aerodynamic forces applied to the wing induce both torsion and bending in the skeleton, but the torsion is greater for the humerus because it is oriented such that forces at the center of pressure have a large torsional moment about the humerus (Swartz et al., 1992⇓). Alteration in the model input parameters can act effectively to shift the position of this center of pressure relative to the humerus, producing the observed sensitivity. The midshaft of the radius, in contrast, is positioned closer to the center of pressure, and small changes in its location have a smaller effect on computed radial stresses.

### Extensions and limitations of the model

There are a number of important limitations to the model in its present
form. It applies strictly only to level, constant-velocity flight of moderate
speed and to a single species, *Pteropus poliocephalus*. With
modifications, this framework can be extended to a greater range of flight
behaviors and to a variety of taxa. However, this will require additional data
and analysis, particularly the detailed documentation of kinematic variation
in relation to activity, body size, wing form, etc., and the detailed analysis
of structural design of the wings of other species. Moreover, as the model is
extended and employed in novel ways, it will not be possible to validate it
*via* the bone strain analysis of level flight in *P. poliocephalus.
In vivo* bone strain measurement is not feasible for very small bones and
so is unlikely to serve as an effective tool to assess flight mechanics in bat
species much smaller than grey-headed flying foxes, i.e. most of the nearly
1000 species of bat. It is conceivable that *in vivo* strain gauge
measurement could be used, however, to obtain data that could serve to
validate models of higher speed or of accelerating or turning flight. It will
be critical, though, to seek additional means of external validation. With
rigorous validation of model accuracy for simple flight behaviors, however, it
is possible to retain a degree of confidence in the model even for behaviors
or taxa for which such detailed validation is impossible.

Finally, we note that this model is not restricted to predicting bone stresses and motions of the animal's center of mass. In future studies, we will use the model to estimate joint forces and moments and flight energetics, as well as extending the model to behaviors including vertical and horizontal acceleration and turning flight. Because the required model inputs, detailed flight kinematics and specific measurements of wing structure and mass distribution are simple in comparison with direct measurement of skeletal stresses, membrane and joint forces, etc., we will be able to extend the model to other taxa. We believe that this kind of model can be employed to provide new insights into flight mechanics that will improve our understanding of wing morphology, kinematics and dynamics. It is a powerful tool for the generation of new hypotheses concerning the mechanics and aerodynamics of bat flight and, by identifying aspects of structural design and flight mechanics that could constrain behavior and influence organismal performance, it may help define and focus future field and morphological studies.

- List of symbols
- A
- area of a single wing
*A*_{c}- cross-sectional area of bone
*A*_{p}- area of a single wing segment
*p* **a**- global acceleration of the shoulder with respect to the non-inertial origin
**a**_{pq}- acceleration on wing segment
*p*at time*q* **a**_{s}- global acceleration of the shoulder relative to the inertial reference frame
*a*_{x.pq}- component of acceleration in the
*x*direction of wing segment*p*at time*q* *a*_{y.pq}- component of acceleration in the
*y*direction of wing segment*p*at time*q*segment*p*at time*q* *a*_{z.pq}- component of acceleration in the
*z*direction of wing segment*p*at time*q* *a*_{sy.q}- vertical acceleration component of the global acceleration at time
*q* *a*_{zy.q}- horizontal acceleration component of the global acceleration at time
*q* - B
- mean wing length
*B*_{mid}- wing length at the middle of the downstroke
*b*_{a}- instantaneous armwing length from shoulder to carpus
*b*_{min}- minimum value of
*b*_{a} - C
- aerodynamic constant
*C*_{D,b}- form drag coefficient of body
*C*_{D,f}- frictional drag coefficient of wing surface
*C*_{m}- added mass coefficient of wing segments
- c
- mean wing chord
*c*_{a}- mean armwing chord
*c*_{p}- wing segment chord
**D**_{b}- body form drag vector
*D*_{b}- magnitude of body form drag
**D**_{f}- skin friction drag vector
*D*_{f}- magnitude of skin friction drag on both sides of both wings
**D**_{i}- mean induced drag vector
*D*_{i}- magnitude of the mean induced drag
**D**_{t}- total drag vector
*D*_{t}- magnitude of the total drag
*d*_{p}- wing segment center of mass position at the middle of downstroke
- E
- complete elliptic integral of the second kind
*E*_{11}- compressional Young's modulus of bone
*E*_{xx}- spanwise plagiopatagium elastic modulus
*E*_{zz}- chordwise plagiopatagium elastic modulus
*e*_{p}- wing segment leading edge position at middle of downstroke
**F**_{pq}- internal force vector on wing segment
*p*at time*q* - F
- magnitude of internal force
*F*_{a}- magnitude of aerodynamic force
**F**_{a.pq}- aerodynamic force vector on wing segment
*p*at time*q* **F**_{d.pq}- net drag vector on wing segment
*p*at time*q* *F*_{d}- magnitude of net drag
**F**_{f}- total force vector acting on skeleton
**F**_{f.pq}- total force vector acting on skeleton due to wing segment
*p*at time*q* *F*_{fx}- magnitude of
*x*component of total force vector**F**_{f} *F*_{fy}- magnitude of
*y*component of total force vector**F**_{f} *F*_{fz}- magnitude of
*z*component of total force vector**F**_{f} *F*_{g}- magnitude of gravitational force
**F**_{g.pq}- gravitational force vector on wing segment
*p*at time*q* **F**_{i.pq}- inertial force vector on wing segment
*p*at time*q* *F*_{i}- magnitude of inertial force
**F**_{l.pq}- lift vector on wing segment
*p*at time*q* *F*_{l}- magnitude of lift
**F**_{m.pq}- added mass force on wing segment
*p*at time*q* *F*_{m}- magnitude of added mass force
**F**_{pq}- internal force
**F**_{s.pq}- total skin internal force on wing segment
*p*at time*q* *F*_{s}- magnitude of total skin internal force on wing segment
*p*at time*q* *F*_{sy.q}- vertical component of the global force at time
*q* **F**_{t.pq}- thrust on wing segment
*p*at time*q* *F*_{t}- magnitude of thrust
*F*_{x}- magnitude of
*x*component of total force vector**F**_{f}in global inertial coordinate system *F*_{y}- magnitude of
*y*component of total force vector**F**_{f}in global inertial coordinate system *F*_{z}- magnitude of
*z*component of total force vector**F**_{f}in global inertial coordinate system *F*_{x′}- magnitude of
*x*component of total force vector**F**_{f}in primed coordinate system *F*_{y′}- magnitude of
*y*component of total force vector**F**_{f}in primed coordinate system *F*_{z′}- magnitude of
*z*component of total force vector**F**_{f}in primed coordinate system - f
- flapping frequency,
*T*^{-1} - g
- gravitational acceleration, 9.81 m s
^{-2} - i
- unit vector of the
*x*axis - J
- polar second moment of area
- j
- unit vector of the
*y*axis - k
- unit vector of the
*z*axis - k
- coefficient of induced drag
- L
- bone length
*L*_{h}- humerus length
*L*_{r}- radius length
- l
- length of wing hairs
*l*_{x}*x*component of moment arm of**F**_{f}at (*x*f,*y*f,*z*f) about (*x*0,*y*0,*z*0)*l*_{y}*y*component of moment arm of**F**_{f}at (*x*f,*y*f,*z*f) about (*x*0,*y*0,*z*0)*l*_{z}*z*component of moment arm of**F**_{f}at (*x*f,*y*f,*z*f) about (*x*0,*y*0,*z*0)- M
- total mass
**M**_{0}- moment vector about the point (
*x*0,*y*0,*z*0) **M**_{0.pq}- moment vector about (
*x*0,*y*0,*z*0) due to forces on wing segment*p*at time*q* *M*_{0x}*x*component of moment of**F**_{f}about (*x*0,*y*0,*z*0)*M*_{0y}*y*component of moment of**F**_{f}about (*x*0,*y*0,*z*0)*M*_{0z}*z*component of moment of**F**_{f}about (*x*0,*y*0,*z*0)*M*_{x′}*x*component of moment of**F**_{f}about (*x*0,*y*0,*z*0) in primed coordinate system*M*_{y′}*y*component of moment of**F**_{f}about (*x*0,*y*0,*z*0) in primed coordinate system*M*_{z′}*z*component of moment of**F**_{f}about (*x*0,*y*0,*z*0) in primed coordinate system*M*_{s}- magnitude of moment of total skin internal force about carpus
- m
- parameter of the complete elliptic integral of the second kind
*m*_{b}- mass of the body alone
*m*_{i}- curve fit coefficient of a wing landmark,
*i*=0-8 *m*_{p}- wing segment mass
- n
- unit vector pointing normal to dorsal face of wing
**n**_{pq}- unit vector pointing normal to dorsal face of wing segment
*p*at time*q* *n*_{x.pq}- component of normal unit vector of wing segment
*p*at time*q*along*x*axis *n*_{y.pq}- component of normal unit vector of wing segment
*p*at time*q*along*y*axis *n*_{z.pq}- component of normal unit vector of wing segment
*p*at time*q*along*z*axis - P
- air pressure
*P*_{c}- proportional location of segmental center of lift relative to chord length
- p
- integer value of wing segment
- Q
- integer time step corresponding to
*t*_{mid} - q
- integer value of time step, ≤40
- R
- radius of bone
*R*_{i}- inner radius of bone
*R*_{ih}- inner radius of humerus
*R*_{ir}- inner radius of radius
*R*_{o}- outer radius of bone
*R*_{oh}- outer radius of humerus
*R*_{or}- outer radius of radius
- S
- mean wing span
- s
- out-of-wing-plane slope of the plagiopatagium at the fifth digit
- T
- wingbeat period
*T*_{b}- magnitude of mean thrust to overcome body form drag
**T**_{b.p}- vector of mean thrust to overcome body form drag on wing segment
*p* *T*_{f}- magnitude of mean thrust to overcome skin friction drag
**T**_{f.p}- vector of mean thrust to overcome skin friction drag on both sides of
wing segment
*p* *T*_{i}- magnitude of the mean thrust to overcome induced drag
**T**_{i.p}- magnitude of the mean thrust to overcome induced drag on wing segment
*p* *T*_{x}- spanwise plagiopatagium tension
*T*_{z}- chordwise plagiopatagium tension
- t
- time, 0≤
*t*≤*T* *t*_{mid}- time at the middle of the downstroke
*T*_{q}- time corresponding to the time step integer
*q* *T*_{s}- plagiopatagium skin thickness
- U
- flight speed
- V
- volume
- W
- shoulder-to-shoulder width of body
*w*_{p}- mean width of a wing segment
*w*_{pq}- width of wing segment
*p*at time*q* - x
- relative coordinate with origin at shoulder, carpus or bone midshaft
*x*_{c}*x*coordinate of the position of the carpus*x*_{f}*x*coordinate of inertial, external or skin membrane force**F**_{f}acting on wing*x*_{0}*x*coordinate of bone location about which moments of**F**_{f}act- x′
- relative coordinate in primed system oriented to bone axes
- Y
- plagiopatagium deflection out of the armwing plane
- y
- relative coordinate with origin at shoulder, carpus or bone midshaft
*y*f*y*coordinate of inertial, external or skin membrane force**F**_{f}acting on wing*y*0*y*coordinate of bone location about which moments of**F**_{f}act- y′
- relative coordinate in primed system oriented to bone axes
- z
- relative coordinate with origin at shoulder, carpus or bone midshaft
*z*f*z*coordinate of inertial, external or skin membrane force**F**_{f}acting on wing*z*0*z*coordinate of bone location about which moments of**F**_{f}act- z′
- relative coordinate in primed system oriented to bone axes
- ϵ1
- longitudinal or normal strain
- ϵ2
- transverse or hoop strain
- ϵxx
- spanwise plagiopatagium strain
- ϵzz
- chordwise plagiopatagium strain
- Γpq
- instantaneous circulation about a wing segment
*p*at time*q* - λ
- angle between armwing and horizontal
*x*axis - νzx
- Poisson's ratio of plagiopatagium due to loading along the
*x*axis - νxz
- Poisson's ratio of plagiopatagium due to loading along the
*z*axis - ρ0
- density of air (1.29 kg m
^{-3}) - σL
- longitudinal or normal stress
- σT
- transverse or hoop stress
- ν1,2
- shear stress
- ψ
- angle between one armwing bone and
*x*′ axis

## ACKNOWLEDGEMENTS

The authors wish to thank Christine Huo for her tireless digitization. We would also like to thank Dr Steven Vogel for allowing P.W. to conduct experiments in the recirculating flume at Duke University. Thanks are also due to Ray and Doralee Watts for logistical support. Dr Geoff Spedding has graciously discussed earlier versions of this work. This work was supported by National Science Foundation grants IBN-9119413 and IBN-9723736 to S.M.S.

- © The Company of Biologists Limited 2001