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

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

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

We employed Newton's second law of motion to balance the inertial force of a wing strip with the external forces as well as the internal force within the wing segment, and then used the segmental balance of forces at each instant in time to solve for the internal force carried by the wing structures. We calculated the stresses developed at the midshaft of the humerus and radius by summing the internal forces over all strips distal to the bone site of interest and computing the moments of these forces about the bone's midshaft. We modified these maximum possible internal forces and moments by subtracting the forces and moments transmitted directly by the plagiopatagium (armwing) to the body, and then employed standard engineering beam theory to convert the remaining three orthogonal forces and moments at the midshaft into normal and shear stresses at the bone surface, and compared stresses calculated for the humeral and radial midshafts with stresses inferred from strains measured empirically in flying bats (Swartz et al., 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.

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

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.

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.

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.

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.

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

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

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

where m₀,..., m₈ 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₀,..., m₅ with which to fit the data. All x direction curve fits have correlation coefficients r≥0.959 whereas all y direction curve fits have correlation coefficients r≥0.993. Displacements in the z direction were an order of magnitude smaller than those in the other two directions and matched the accuracy of the position data in all three directions, resulting in curve fits with correlation coefficients r>0.44.

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

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

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.

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_pqc_p²)/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_pqc_p²ρ0)/4 or m_p≈0.6 g for P. poliocephalus.

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.

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

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

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

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

where W is the shoulder-to-shoulder width of the body, and πW²/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 Mg is total weight, k=1.2 is a constant proposed by Pennycuick (Pennycuick, 1989) and S is the mean wing span. Induced drag acts in the opposite direction to the direction of flight. If the instantaneous wing span is defined as the shortest distance from one wingtip to the other along the wing surface, then one standard deviation in the wing span in 7.5% of the mean.

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.

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

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

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:

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

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.

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

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

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

where W is body width and F_a is the total aerodynamic force, a simple function of total drag and body weight.

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

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:

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.

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.

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

Although inertial forces exerted on the skin can probably approach the magnitude of aerodynamic forces near the carpus during periodic rapid wing acceleration, we modeled only the skin deformation induced by the instantaneous lift acting on the plagiopatagium. We divided total armwing lift by the instantaneous plagiopatagium area b_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:

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²/(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_xc_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.5c_a in this special coordinate system.

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

Assuming that a bending moment induces a linear increase in normal stress σ _L with increasing distance from the neutral axis, we calculated the magnitude of the normal stress induced by M_y′ as:

where J, the polar second moment of area of the bone cross section, is π(R_o⁴-R_i⁴), and we are free to choose R_i≤z′≤R_o. For positive M_y′, the cranial face of the bone is in tension while the caudal face is in compression, as indicated by the sign of σ _L. The magnitude of the normal stress induced by M_z′ is:

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

where R_i≤R≤R_o. This shear stress is uniform around the bone circumference, and the maximum values of both the normal and shear stresses occur at the bone outer radius R=R_o.

The normal stress induced by the force F_x was calculated by assuming that the force is distributed evenly over the cross section, giving rise to a normal stress:

where A_c is the bone's cross-sectional area. Forces acting perpendicular to the long axis of the bone give rise to shear stresses that are not distributed evenly along the bone cross section (Beer and Johnston, 1981). The magnitude of the shear stress at y′=0, z′=±R_o due to the internal force F_y′ is:

where, for positive F_y′, shear stress is positive on the cranial (leading) and negative on the caudal (trailing) surface of the bone. The magnitude of the shear stress due to the internal force F_z′ at z′=0, y′=±R_o is:

where, for positive F_z′, shear stress is positive on the dorsal and negative on the ventral surface of the bone. The two magnitudes of τ _1,2 given here apply to only two points, each along the outer perimeter of the bone.

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

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

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

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

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

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

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

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.

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

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

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

Empirically recorded bone strains indicated relatively large transverse stresses at the humeral and radial midshafts (Swartz et al., 1992EF57). The constitutive model for bone presented by Carter (Carter, 1978EF11) is almost volume-conserving and predicts ϵ₂=0.46ϵ₁ 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.

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

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.

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

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

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

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

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

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.

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

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

List of symbols
 
• A

  area of a single wing

 
• A_c

  cross-sectional area of bone

 
• A_p

  area of a single wing segment p

 
• a

  global acceleration of the shoulder with respect to the non-inertial origin

 
• a_pq

  acceleration on wing segment p at time q

 
• a_s

  global acceleration of the shoulder relative to the inertial reference frame

 
• a_x.pq

  component of acceleration in the x direction of wing segment p at time q

 
• a_y.pq

  component of acceleration in the y direction of wing segment p at time q segment p at time q

 
• a_z.pq

  component of acceleration in the z direction of wing segment p at time q

 
• a_sy.q

  vertical acceleration component of the global acceleration at time q

 
• a_zy.q

  horizontal acceleration component of the global acceleration at time q

 
• B

  mean wing length

 
• B_mid

  wing length at the middle of the downstroke

 
• b_a

  instantaneous armwing length from shoulder to carpus

 
• b_min

  minimum value of b_a

 
• C

  aerodynamic constant

 
• C_D,b

  form drag coefficient of body

 
• C_D,f

  frictional drag coefficient of wing surface

 
• C_m

  added mass coefficient of wing segments

 
• c

  mean wing chord

 
• c_a

  mean armwing chord

 
• c_p

  wing segment chord

 
• D_b

  body form drag vector

 
• D_b

  magnitude of body form drag

 
• D_f

  skin friction drag vector

 
• D_f

  magnitude of skin friction drag on both sides of both wings

 
• D_i

  mean induced drag vector

 
• D_i

  magnitude of the mean induced drag

 
• D_t

  total drag vector

 
• D_t

  magnitude of the total drag

 
• d_p

  wing segment center of mass position at the middle of downstroke

 
• E

  complete elliptic integral of the second kind

 
• E₁₁

  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 (xf, yf, zf) about (x0, y0, z0)

 
• l_y

  y component of moment arm of F_f at (xf, yf, zf) about (x0, y0, z0)

 
• l_z

  z component of moment arm of F_f at (xf, yf, zf) about (x0, y0, z0)

 
• M

  total mass

 
• M₀

  moment vector about the point (x0, y0, z0)

 
• M_0.pq

  moment vector about (x0, y0, z0) due to forces on wing segment p at time q

 
• M_0x

  x component of moment of F_f about (x0, y0, z0)

 
• M_0y

  y component of moment of F_f about (x0, y0, z0)

 
• M_0z

  z component of moment of F_f about (x0, y0, z0)

 
• M_x′

  x component of moment of F_f about (x0, y0, z0) in primed coordinate system

 
• M_y′

  y component of moment of F_f about (x0, y0, z0) in primed coordinate system

 
• M_z′

  z component of moment of F_f about (x0, y0, z0) 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₀

  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

 
• yf

  ycoordinate of inertial, external or skin membrane force F_f acting on wing

 
• y0

  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

 
• zf

  z coordinate of inertial, external or skin membrane force F_f acting on wing

 
• z0

  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 thex axis

 
• νxz

  Poisson's ratio of plagiopatagium due to loading along thez 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

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.