QUICK SEARCH:   [advanced] Author: Keyword(s):
Home     Help     Feedback     Subscriptions     Archive     Search     Table of Contents

This Article Summary Figures Only Full Text (PDF) Alert me when this article is cited Alert me if a correction is posted Services Email this article to a friend Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via Google Scholar Google Scholar Articles by WATTS, P. Articles by SWARTZ, S. M. Search for Related Content PubMed PubMed Citation Articles by WATTS, P. Articles by SWARTZ, S. M.

The Journal of Experimental Biology 204, 2873-2898 (2001)
© 2001 The Company of Biologists Limited

A COMPUTATIONAL MODEL FOR ESTIMATING THE MECHANICS OF HORIZONTAL FLAPPING FLIGHT IN BATS : MODEL DESCRIPTION AND VALIDATION

PHILIP WATTS¹, ERIKA J. MITCHELL¹ and SHARON M. SWARTZ^1,2,*

¹ Department of Ecology and Evolutionary Biology, Brown University, Providence, RI 02912, USA
² Division of Engineering, Brown University, Providence, RI 02912, USA
^* Author for correspondence (e-mail: sharon_swartz{at}brown.edu )

Accepted June 5, 2001

	Summary

TOP Summary Introduction Materials and methods Results Discussion References

 
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.

Key words: flight, wing mechanics, aerodynamics, computational modelling, Chiroptera, bat, Pteropus poliocephalus

	Introduction

TOP Summary Introduction Materials and methods Results Discussion References

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

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

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

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

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

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

	Materials and methods

TOP Summary Introduction Materials and methods Results Discussion References

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

View larger version (33K): [in this window] [in a new window]   Fig. 1. Plan-view of the ventral surface of the wing of a Pteropus poliocephalus held in a horizontal plane, a position similar to that at the middle of the downstroke. The wing is subdivided into 14 chordwise strips; there are seven segments between the shoulder and the carpus (the plagiopatagium or armwing), three between the carpus and the metacarpophalangeal (MCP III) joint of the third digit (the proximal handwing) and four between the MCP III joint and the wingtip (the distal handwing). The locations of the centers of mass of these strips relative to a reference line connecting the two shoulders (broken line) are indicated by the filled circles. The large circles labeled shoulder (glenohumeral joint), carpus, MCP III and wingtip were used as digitizing markers for collecting kinematic data. For the third wing segment (subscript p), the length of the wing chord, cp, the leading edge position, ep, and the distance from the segmental center of mass to the reference line, dp, are also indicated. The variation in these parameters within the wingbeat cycle is represented in the model by a second subscript, q, that ranges from 1 to 40 with the 40 equal time increments within the complete wingbeat cycle.

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

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

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

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

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

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

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

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

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

View larger version (13K): [in this window] [in a new window]   Fig. 2. Orientations of and relationships between the inertial and non-inertial coordinate systems; the actual and mean paths of the shoulder are depicted from a lateral view of the bat moving from left to right across the figure. The x axis is perpendicular to the plane of the page and is represented by the filled circles.

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

View larger version (16K): [in this window] [in a new window]   Fig. 3. Representative illustration of the complex kinematics of the bat wing after curve-fitting raw data; this view depicts the movements of the carpus (open squares) and wingtip (filled circles) with respect to the shoulder joint during both a downstroke and an upstroke. The lines connecting the shoulder (coordinates 0,0) and the carpus schematically represent armwing length, and those between the carpus and the wingtip, handwing length. Time intervals between data points are equal and show that the wrist, in particular, moves far more rapidly in midstroke than at other times in the wingbeat cycle.

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:

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

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

Gravitational force
The action of gravity on the wing of a large bat is typically an order of magnitude smaller than the aerodynamic force encountered in level flight because the former acts only on the wing mass while the latter supports the entire body mass (e.g. Thollesson and Norberg, 1991). Nevertheless, the segmental gravitational force is given by:

(6)

where m_p is the segmental mass and j is a unit vector of the global y axis pointing in the direction opposite to gravitational acceleration, g. The gravitational force acts at the center of mass of each wing segment (see Table 2 for the chordwise positions of the centers of mass relative to the positions of the leading edge). We assumed that the center of mass lies at midspan of each wing segment.

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

Inertial force
Inertial forces on the wing vary with wing motion and thus can range in magnitude from zero during episodes of gliding to values comparable with the large aerodynamic forces exerted on the wing during flapping flight. By definition, the inertial force also represents the net sum of all external forces as well as the force carried by the wing structures for a given wing segment. On the basis of segmental accelerations, the inertial force on wing segment p at time q is given by:

(7)

where each component of a force is treated independently by the model. The inertial force acts at the center of mass and in the same vector direction as the acceleration of a wing segment. The horizontal acceleration s_sz.q and the vertical acceleration s_sy.q are the two components of the global acceleration at each instant of time and must be solved for iteratively in the program. We note that the Einstein summation convention is not implied in this paper by the use of repeated subscripts.

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

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

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

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,b0.5 (White, 1991). In contrast to the bodies of birds, those of bats are not streamlined, and can, indeed, resemble a sphere to some extent because of the proportions of the ribcage. Moreover, form drag is smaller than the induced drag, which mitigates against searching for a precise numerical value of the drag coefficient. Our results for this species (see below) demonstrate that there is not much evolutionary pressure to drive bat body streamlining from a biomechanical perspective. It follows that the form drag D_b is:

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

The total drag D_t on the bat in forward flight is the sum of skin friction drag, body form drag and induced drag:

(13)

where D_t is the magnitude of the total drag. We have assumed that the total drag is constant over time because it is an order of magnitude smaller than the aerodynamic force. On the basis of the data in Table 1, the ratio of mean total lift Mg to total drag is Mg/D_t9 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_szD_t/M over a duration T/4, then we can write the fractional change in flight velocity as approximately:

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

(16)

(17)

to be consistent with the value of skin friction drag above. We calculate that A=0.0633 m² 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, 2B. We write the mean thrust needed to overcome form drag and induced drag on a wing segment p as:

(20)

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

Passive aeroelastic wings supported by leading edge spars tend to produce maximum thrust during the middle of the downstroke as a result of wing twisting that rotates the aerodynamic force towards the direction of flight (DeLaurier, 1993c). Therefore, we assumed sinusoidal variation in thrust over the wingbeat cycle, with the maximum thrust occurring during the middle of the downstroke. We calculated the magnitude of the thrust explicitly to find the orientation of the aerodynamic force. In the absence of specific thrust observations, we verified a posteriori that sinusoidal thrust variation does not strongly affect our simulation results. Ideally, a more precise thrusting function is desirable, particularly one that captures the 65% downstroke versus 35% upstroke cycle, but we did not feel that such an effort was germane to this work at present. We write the magnitude of the thrust at wing segment p and time q as:

(21)

where f=1/T is the flapping frequency, the argument within the cosine function is expressed in radians, and the time t_q associated with the integer q is written:

(22)

where t_mid is the time at the middle of the downstroke, and t=T/40 is the time step. The maximum thrust 2(T_f.p+T_b.p+T_i.p) occurs when tqt_mid, and thrust is equal to zero when tqt_mid+(T/2) during the upstroke. The mean thrust generated by both wings over one wingbeat cycle is equal to the total drag since the contribution from the cosine function averages to zero.

Aerodynamic force distribution and constant-circulation flight
The magnitude of the total aerodynamic force F_a needed to sustain forward flight is a combination of weight and total drag components with a mean value over one wingbeat cycle of:

(23)

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

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

In our model, we defined the magnitude of the aerodynamic force on a wing segment p as:

(24)

where _pq is the circulation about the midline of the wing segment and w_pq is the instantaneous width of the wing segment. We invoked the constant-circulation hypothesis when we made the magnitude of the aerodynamic force _pqw_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_qt_mid=00.38T 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:

View larger version (20K): [in this window] [in a new window]   Fig. 4. Elliptical distribution of the magnitude of aerodynamic forces at each wing segment; segment 1 is adjacent to the shoulder and segment 14 includes the wingtip (see Fig. 1).

Net drag and lift
We defined net drag by adding the drag and thrust acting on each wing segment. The net drag on wing segment p at time q is given by:

(28)

and acts in the plane of the wing segment opposing the direction of flight. We calculated the magnitude of the lift from the mathematical relationship between the magnitude of the aerodynamic force and the magnitude of its two orthogonal components, lift and thrust. Hence, the lift on wing segment p at time q is written:

Internal force
The internal force is the force carried by the wing structures that enables the wing to undergo the observed accelerations and to resist the external forces applied to the wing. We invoke the common form of Newton's second law of motion for wing segment p at time q so that:

(30)

where the internal force F_pq is not given any subscript. We did not calculate the internal force explicitly in the model since we were interested in the relative contribution of inertial and external forces to bone stresses. Moreover, the inertial and external forces act at different locations on a given wing segment, requiring both a net internal force and a net internal moment about the center of mass for each wing segment. Instead, we relied on the definition of the internal force as the inertial force minus the gravitational force, added mass force, lift and net drag, and computed the contributions from these forces separately. We also calculated the component of the internal force associated with tension in the plagiopatagium (see below, section on plagiopatagium tension).

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

(34)

where a_sy.q is the vertical component of the global acceleration, and F_a.0 is the aerodynamic force generated by the body. The horizontal force balance provides at time q:

(35)

where a_sz.q is the horizontal component of the global acceleration. Form and induced drag of the body have already been accounted for in the net drag of each wing segment. We solved these two equations for the two components of the global acceleration iteratively because the equations for F_y.pq and F_z.pq contain the global acceleration components within inertial terms.

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

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^-1s 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.

View larger version (17K): [in this window] [in a new window]   Fig. 5. Schematic diagram of the rectangular idealization of the armwing membrane. The x and z coordinate axes used to describe membrane deflection mathematically are indicated; the y axis is perpendicular to the plane of the membrane. Fs indicates the point of application and approximate orientation of the skin internal force, ca is the mean armwing chord and ba is the instantaneous shoulder-to-carpus distance.

We estimated the spanwise (E_xx) and chordwise (E_zz) elastic moduli of P. poliocephalus plagiopatagia from published mean values of moduli of plagiopatagia of other bats, yielding E_xx =3 MPa and E_xzz=37 MPa; these values correspond closely to the elastic moduli found for the plagiopatagia of A. jamaicensis, the only large fruit-eating bat from which these data have been collected (Swartz et al., 1996). However, because this species is an order of magnitude smaller in body mass than P. poliocephalus, we also explored the effects of lower moduli on force transfer. We further assumed the Poisson's ratio _zx=1 as a typical value for skin composed of crossed fiber sheaths (Frolich et al., 1994). Symmetry of the skin stress tensor associated with mechanical equilibrium requires that _xz=(E_xxzx)/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_ac_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:

(38)

and

(39)

where we assumed uniform plane stress across the membrane thickness t_s=0.2 mm. We evaluated the spanwise (_xx) and chordwise (_zz) skin strains by assuming that, when under no tension, the membrane had a chordwise length c_a and a spanwise length b_min, the minimum value of b_a. It follows from the definition of an arc length that the strain _xx is:

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

Forces, moments and stresses on the wing skeleton
To analyze the effect of flight-related forces on the skeleton, we designated the point at which an inertial or external or skin membrane force vector F_f acts on a given wing segment as (x_f, y_f, z_f) and the point on the skeleton about which we seek the moment M₀ as (x₀, y₀, z₀). By varying the location of the point for which the analysis is carried out, we may then analyze the skeletal loading and bending at any location on the wing skeleton of the armwing that is biologically significant and through which a known fraction (usually taken as 100%) of the remaining internal force is transmitted. Although the analysis applies generally to any location along the armwing, we focused here on the moments about the humeral and the radial midshafts, to carry out comparisons of model estimates with empirically measured values (Swartz et al., 1992). We defined three moment arms: