## 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., 1989EF27; Swartz, 1997EF55; Swartz, 1998EF56; Thomas, 1975EF61; Thomas and Suthers, 1972EF62; Winter et al., 1993EF67), 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, 1995EF49), 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., 1972EF22; Hermanson and Alternbach, 1983EF24; Holbrook and Odland, 1978EF25; Norberg, 1970aEF30; Norberg, 1972aEF32; Norberg, 1972bEF33; Papadimitriou et al., 1996EF38; Strickler, 1978EF53; Swartz, 1997EF55; Swartz, 1998EF56; Swartz et al., 1992EF57; Swartz et al., 1996EF59; Vaughan, 1959EF63; Vaughan, 1970bEF65). Similarly, significant contributions have been made to our understanding of the kinematics of bat flight (e.g. Aldridge, 1986EF2; Aldridge, 1987EF3; Baggøe, 1987EF4; Brandon, 1979EF9; Norberg, 1970bEF31; Norberg, 1976aEF34; Norberg 1976bEF35; Norberg, 1990EF36; Rayner, 1987EF46; Rayner et al., 1986EF47; Vaughan, 1970aEF64). 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, 1972aEF32; Pennycuick, 1973EF40), the velocity-dependence of wing kinematics (Aldridge, 1986EF2) and the relationship between wing inertia, energetics and maneuverability (Norberg, 1976aEF34; Thollesson and Norberg, 1991EF60).

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, 1977EF7; DeLaurier, 1993aEF14; Ellington, 1975EF18; Ellington, 1978EF19; Ellington, 1984EF20; Ennos, 1988EF21; Norberg, 1975EF29; Norberg, 1970bEF31; Norberg, 1976aEF34; Pennycuick, 1968EF39; Pennycuick, 1975EF41; Rayner, 1986EF45; Rees, 1975EF48; Spedding, 1992EF51; Spedding and DeLaurier, 1995EF52; Withers, 1981EF68).

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., 1992EF57). 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, 1973EF5).

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, 1973EF5) 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., 1992EF57). 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, 1985EF10; Swartz, 1998EF56; Swartz et al., 1992EF57; Swartz et al., 1993EF58). 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, 1973EF5). 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, 1991EF17).

### 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., 1992EF57; Swartz et al., 1993EF58). 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, 1978EF11; Beer and Johnston, 1981EF6; Biewener, 1983EF8; Currey, 1987EF13). 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, 1993bEF15; Norberg, 1976aEF34; Thollesson and Norberg, 1991EF60). 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, 1991EF60) (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, 1973EF5). 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., 1992EF57; Swartz et al., 1993EF58). 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., 1989EF44). 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

*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,

**. 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.**

*g*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

*p*at time

*q*is given by:

*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, 1990EF36) at the middle of the downstroke as a plausible approximation for a bat wing in fast forward flight (Rayner et al., 1986EF47; Rayner, 1986EF45; Spedding, 1987EF50) (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, 1995EF52); 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, 1993bEF15); 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, 1972aEF32), 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, 1991EF66). 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, 1993bEF15). 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, 1990EF36). The induced drag **D**_{i} becomes:

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

*k*=1.2 is a constant proposed by Pennycuick (Pennycuick, 1989EF42) 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.

**D**

_{t}on the bat in forward flight is the sum of skin friction drag, body form drag and induced drag:

*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*

**to total drag is**

*g**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, 1995EF52). 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:

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:

*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.

*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:

*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:

*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:

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, 1993bEF15). 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.

*p*as:

_{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, 1990EF36), 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

*p*at time

*q*is given by:

*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, 1991EF26).

### Internal force

*p*at time

*q*so that:

**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:

*q*is:

*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*:

*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., 1996EF59). 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., 1994EF23). 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:

*T*

_{x}and

*T*

_{z}are functions of

*Y*through the orthotropic constitutive relationships:

*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

**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., 1992EF57). We defined three moment arms:

**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:

**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, 1991EF60).

*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:

*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., 1992EF51).

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, 1981EF6). 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., 1992EF57) assumed that compact cortical bone is a transversely orthotropic material and used published values for compact bone moduli and Poisson's ratio (Carter, 1978EF11) 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, 1984EF12). 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, 1972EF54; Swartz et al., 1996EF59). 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., 1996EF59), 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., 1992EF57), 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., 1992EF57). The constitutive model for bone presented by Carter (Carter, 1978EF11) 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., 1992EF57). 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, 1972EF54; Swartz et al., 1996EF59), 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., 1996EF59), 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., 1992EF57). 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.

- 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.

## References

**Abramowitz, M. and Stegun, I. A.**(

**Aldridge, H. D. J. N.**(

*Rhinolophus ferrumequinum*, in horizontal flight at various flight speeds.

**Aldridge, H. D. J. N.**(

**Baggøe, H. J.**(

**Bartholomew, G. A. and Carpenter, R. E.**(

**Beer, F. P. and Johnston, E. R., Jr**(

**Bennett, L.**(

**Biewener, A. A.**(

**Brandon, C.**(

**Carpenter, R. E.**(

*Pteropus poliocephalus.*

**Carter, D. R.**(

**Currey, J. D.**(

**Currey, J. D.**(

**DeLaurier, J. D.**(

**October**,

**DeLaurier, J. D.**(

**April**,

**DeLaurier, J. D.**(

**May**,

**Eby, P.**(

*Pteropus poliocephalus*(Chiroptera: Pteropodidae), from two maternity camps in northern New South Wales.

**Ellington, C. P.**(

*Encarsia formosa*. In

**Ellington, C. P.**(

**Ellington, C. P.**(

**Ennos, A. R.**(

**Findley, J. S., Studier, E. H. and Wilson, D. E.**(

**Frolich, L. M., LaBarbera, M. and Stevens, W. P.**(

**Hermanson, J. W. and Alternbach, J. S.**(

*Artibeus jamaicensis.*

**Holbrook, K. and Odland, G. F.**(

**Katz, J. and Plotkin, A.**(

**Kurta, A., Bell, G. P., Nagy, K. A. and Kunz, T. H.**(

*Myotis lucifugus*).

**Meriam, J. L. and Kraige, L. G.**(

**Norberg, R. A.**(

*Aeschna juncea*L., kinematics and aerodynamics. In

**Norberg, U. M.**(

*Plecotus auritus*Linnaeus (Chiroptera).

**Norberg, U. M.**(

*Plecotus auritus*Linnaeus.

**Norberg, U. M.**(

**Norberg, U. M.**(

*Rousettus aegypticus*(E. Geoffroy) (Pteropidae).

**Norberg, U. M.**(

*Plecotus auritus.*

**Norberg, U. M.**(

*Plecotus auritus.*

**Norberg, U. M.**(

**Norberg, U. M. and Rayner, J. M. V.**(

**Papadimitriou, H. M., Swartz, S. M. and Kunz, T. H.**(

*Tadarida brasiliensis.*

**Pennycuick, C. J.**(

*Columba livia.*

**Pennycuick, C. J.**(

**Pennycuick, C. J.**(

**Pennycuick, C. J.**(

**Popov, E. P.**(

**Press, W. H., Flannery, B. P., Teukolsky, S. A. and Vetterling, W. T.**(

**Rayner, J. M. V.**(

**Rayner, J. M. V.**(

**Rayner, J. M. V., Jones, G. and Thomas, A.**(

**Rees, C. J. C.**(

**Simmons, N. B.**(

**Spedding, G. R.**(

*Falco tinnunculus*) in flapping flight.

**Spedding, G. R.**(

**Spedding, G. R. and DeLaurier, J. D.**(

**Strickler, T. L.**(

**Studier, E. H.**(

**Swartz, S. M.**(

**Swartz, S. M.**(

**Swartz, S. M., Bennett, M. B. and Carrier, D. R.**(

**Swartz, S. M., Bennett, M. B., Gray, J. A. and Parker, A.**(

*in vivo*loading in the distal wing of fruit bats.

**Swartz, S. M., Groves, M. D., Kim, H. D. and Walsh, W. R.**(

**Thollesson, M. and Norberg, U. M.**(

**Thomas, S. P.**(

*Phyllostomus hastatus*and

*Pteropus gouldii.*

**Thomas, S. P. and Suthers, R. A.**(

**Vaughan, T. A.**(

**Vaughan, T. A.**(

**Vaughan, T. A.**(

**White, F. M.**(

**Winter, Y., Helverson, O. V., Norberg, U., Kunz, T. H. and Steffensen, J. F.**(

*Glossophaga soricina*(Phyllostomidae: Glossophaginae). In

**Withers, P. C.**(