Flexural stiffness in insect wings II. Spatial distribution and dynamic wing bending

During flapping flight, insect wings accelerate masses of air, generating the forces necessary to support the insect's weight and to perform complex maneuvers. At the same time, these ultralight airfoils (generally only 0.5-5%of body mass; Ellington,1984b) must withstand the forces imposed upon them by the surrounding air, as well as the inertial forces caused by accelerating and decelerating their own mass up to several hundred times per second.

Insect wings perform these roles extremely successfully, despite the fact that they are largely passive structures, with no muscular control past the wing base (Wootton, 1992). Although they are strengthened by a network of tubular veins, the wings of many species deform noticeably during flight, especially during slow flight and hovering (Willmott and Ellington,1997a). These dynamic changes in the three-dimensional shape of wings could potentially affect many aspects of force production, yet few models of insect flight have successfully incorporated passive wing flexibility or examined the effects of flexibility on force production.

In the past two decades, models of insect flight, both mathematical and physical, have contributed enormously to our understanding of the basic mechanisms of force production, despite assuming that insect wings are rigid structures [e.g. mathematical models of Ellington(1984a), Lighthill(1973), Sane and Dickinson(2002), Savage et al. (1978),Smith et al. (1996) and Wilkin and Williams (1993), and physical models of Bennett(1966, 1970), Dickinson et al.(1999), Ellington et al.(1996) and Spedding and Maxworthy (1985)]. Some aspects of wing flexibility have been mimicked in models by altering the relative positions of wing regions(Liu et al., 1998; Vest and Katz, 1996) or by modeling deformations as harmonic waves(Combes and Daniel, 2001; Daniel, 1987; Wu, 1971). These approaches provide unique insights into the mechanisms of force generation during flight,but often neglect one or more critical components of wing deflection (e.g. spanwise bending, chordwise bending or torsion), which can have large effects on aerodynamic force production (Batchelor,1967). Models of insect flight that incorporate passive wing flexibility (in which shape changes are driven by forces imposed upon the wing rather than being specified in advance) are exceedingly rare (e.g. Smith, 1996).

One difficulty in modeling passive wing flexibility is that forces applied at the wing base lead to bending and twisting that are influenced not only by overall stiffness and gross anatomical features (e.g. flexion or fold lines),but also by the distribution of flexibility throughout the wing(Wootton, 1999). Extremely detailed finite element models of insect wings (incorporating vein configuration, three-dimensional relief, and variations in vein and membrane characteristics) could potentially reproduce the distribution of flexibility in real wings by accounting for the precise structural and material properties of the wing. Unfortunately, precise data about the local properties of insect wings are often unavailable and finite element models must be constructed using simplifying assumptions. For example, many finite element models assume that the material stiffness of the wing (Young's modulus, E) is equivalent to that of insect chitin and constant throughout the wing(Smith, 1996; Kesel et al., 1998); however,recent measurements reveal that E can vary widely within a wing(Smith et al., 2000) and that other proteins, such as resilin, occupy key positions (e.g. wing vein joints)in insect wings (Gorb, 1999; Haas et al.,2000a,b). Measuring spatial variation in these material properties, as well as the details of vein and membrane structure, is a time-consuming process that would need to be repeated for each new species studied.

Although these detailed approaches can provide important information about functional wing morphology in pivotal, well-studied species (e.g. Herbert et al., 2000), a more general approach to wing stiffness measurements and modeling could facilitate comparative studies and attempts to incorporate passive flexibility into models of force production. Rather than measuring (and modeling) geometric and material properties separately, the overall bending response of the wing, or flexural stiffness (EI, the product of material stiffness Eand second moment of area I) can be determined. Performing measurements of EI averaged over the whole wing is relatively straightforward (see Combes and Daniel,2003a), but measuring spatial variation in flexural stiffness throughout a wing is more challenging. Steppan(2000) approached this issue by measuring average flexural stiffness over increasingly larger sections of dried butterfly wings, and Wootton et al.(2000) measured average flexural stiffness in three isolated sections of locust hindwing. However, no studies to date have demonstrated a method of measuring the spatial variation in local flexural stiffness of intact insect wings.

In this study, we developed a method of approximating spatial variation in flexural stiffness along two axes of the wing (in the spanwise and chordwise direction). We measured flexural stiffness variation in two insects, the hawkmoth Manduca sexta and the dragonfly Aeshna multicolor. These insects have wings of similar size and are both agile fliers, capable of hovering as well as fast, forward flight. However, they are very distantly related and display large differences in wing shape and venation pattern (see Combes and Daniel, 2003a) that may underlie differences in the distribution of flexural stiffness in their wings.

To determine the spatial pattern of flexural stiffness in wings, we first developed a method to measure displacement (due to a point force) along the wing in the spanwise and chordwise directions. Next, we proposed various alternatives for how flexural stiffness might vary along the wing (represented by simple mathematical functions) and predicted the patterns of displacement that a loaded wing with these theoretical stiffness distributions would display. We then found the flexural stiffness distribution that produced a pattern of displacement most similar to the pattern measured in a real wing. We used this method to estimate the flexural stiffness distribution of the wing in response to forces applied from both the dorsal and the ventral side,as dorsal/ventral asymmetry has been noted in previous studies(Ennos, 1988; Steppan, 2000; Wootton, 1993; Wootton et al., 2000).

We also created two different finite element models of an insect wing,based on the planform geometry of a Manduca sexta forewing. Rather than realistically reproducing the structure and behavior of Manducawings, these generalized models provide a method of assessing the consequences of flexural stiffness distributions to wing bending. We attempted to capture major elements of wing structure in the models (e.g. vein position and flexural stiffness distribution) while avoiding detailed structural features,such as variation in wing thickness and three-dimensional relief. We used these models to examine how different spatial patterns of EIdetermine bending patterns in wings subjected to both static and dynamic loads.

Insect handling and wing preparation were performed as in Combes and Daniel(2003a). Calculations of EI distribution require both the applied force and the displacement along the length of the wing. We measured displacement of wings in the spanwise direction along a line running from the wing base to the tip, and in the chordwise direction along a line running from the leading edge to the trailing edge (approximately midway between the base and tip). We attached a fresh wing at either the wing base or leading edge, and applied a point force to the wing with a pin attached to a force beam (see Combes and Daniel, 2003a),using a small drop of cyanoacrylate glue to prevent the pin from slipping off the wing tip or trailing edge. We illuminated the wing with sheets of laser light; these sheets appear as lines on the wing, and can be adjusted so that one line illuminates a strip of the wing from the base to the tip(Fig. 1A) or from the leading to the trailing edge (intersecting the pin). When a force is applied to the wing, the laser lines follow the surface of the wing, shifting their position in proportion to the displacement of the wing(Fig. 1B). We used a Coolpix 900 digital camera (Nikon, Tokyo, Japan) angled approximately 45° from the surface of the wing to photograph unloaded and loaded wings, as well as a calibration object containing a series of steps of known height.

We measured displacement in response to loads of varying magnitude,returning to the unloaded position to photograph the wing before applying each new load. We checked the repeatability of displacement measurements on the wings of a hawkmoth (in both the spanwise and chordwise direction) by applying five different loads to the dorsal surface of the wing and repeating each load three times (applying loads in random order).

For all other hawkmoth and dragonfly wings, we measured displacement on the dorsal surface of the wing at four different loads, gently removed the glued pin from the edge of the wing (without tearing the wing membrane), flipped the slide over, and measured displacement on the ventral side of the wing at four different loads. Dorsal stiffness in this study describes the stiffness of the wing in response to loading on the dorsal surface (resulting in a dorsally convex surface; see Steppan,2000), which is equivalent to ventral flexion in other studies(Ennos and Wootton, 1989; Wootton, 1981, 1992).

After completing the flexural stiffness measurements, we photographed the wings on a white background from above and measured wing area, span and maximum chord length in NIH Image. For Manduca sexta (whose wing morphology may depend on gender; Willmott and Ellington, 1997b), we also calculated aspect ratio(span²/area) of the forewings and wing loading (body mass divided by the area of both forewings).

In each hawkmoth, we measured spanwise flexural stiffness on one forewing and chordwise flexural stiffness on the other. We did not measure the flexural stiffness of the hindwings, which are much smaller and overlap with the forewings to varying degrees during flight. In dragonflies (which have independent fore- and hindwings of similar size), we measured spanwise flexural stiffness on one forewing and one hindwing, and chordwise flexural stiffness on the other forewing and hindwing.

We analyzed images of loaded and unloaded wings with a custom Matlab program (developed by A. Trimble) that finds the center of the laser line running from the wing base to the tip or from the leading to the trailing edge(Fig. 1A,B). We used a 2^nd order Butterworth filter to remove noise from the line position data, and splined the unloaded and loaded data to an equal number of points for comparison. We found the difference in line position between the two data sets [y(x); Fig. 1B] and converted this difference to actual wing displacement along the wing δ(x) with a factor derived from analysis of the calibration object.

With this displacement data and the applied force, we were able to estimate local flexural stiffness by using a continuous beam equation to approximate EI variation along the wing (see Appendix). This equation defines local flexural stiffness as a function of the local curvature (second spatial derivative of displacement) and the local moment (applied force times distance from the point of force application). To avoid errors caused by differentiating displacement data, we solved this as an inverse problem; we posed several simple mathematical functions that might approximate flexural stiffness variation along the wing (constant, linear, exponential, or 2^nd degree polynomial) and calculated the expected pattern of displacement due to the applied force and these possible EIdistributions. We then used a simplex minimization program to find the EI distribution that provided the best fit to the measured displacement along the wing. Finally, we calculated average flexural stiffness(⁠

To verify the method of determining flexural stiffness distribution, we measured δ(x) and calculated local flexural stiffness of a rectangular glass coverslip, which is made of homogeneous material and has a constant material stiffness (E) and second moment of area(I).

Because wing displacement in hawkmoth and dragonfly wings does not increase linearly with applied force (see Combes,2002), our estimates of flexural stiffness varied with the relative displacement of the wing (δ_T/L, tip displacement divided by wing span or trailing edge displacement divided by chord length). To facilitate comparisons between individuals and species, we standardized measurements at a relative displacement of 0.05 for spanwise measurements and 0.08 for chordwise measurements (values that are within the range measured on real wings). For each individual, we used the slope of the relationship between EI variables (⁠

In Manduca sexta, we collected spanwise data from nine males and four females, and chordwise data from ten males and nine females. In Aeshna multicolor, we used only males, but tested both fore and hindwings (which have different morphologies), so we examined stiffness separately in the two sets of wings. We collected spanwise and chordwise data on both fore and hindwings from eight individuals.

Within each group, we tested for differences between the dorsal and ventral side of the wing with a Wilcoxon signed rank test. We also tested for differences between male and female Manduca sexta with a Mann-Whitney U-test.

To investigate how flexural stiffness variation affects wing bending, we created two simplified finite element models of an insect wing (based on the forewing of a male Manduca sexta) with the same geometry, but dramatically different patterns of spatial variation in flexural stiffness(Fig. 2). We used MSC Marc/Mentat 2001 to construct model wings with an accurate planform configuration of veins and membrane, but a simplified three-dimensional geometry (with no camber or three-dimensional relief; see Discussion). We increased the material stiffness of vein elements beyond that of the surrounding membrane to mimic the higher second moment of area of tubular veins and produce the measured spanwise-chordwise anisotropy in flexural stiffness (see Combes and Daniel,2003a).

The models were composed of thin shell elements with a density of 1200 kg m^-3 (as measured in insect wings; Wainwright et al., 1982) and a thickness of 45 μm. We used a Poisson's ratio of 0.49, as measured in some biological materials (Wainwright et al.,1982); because the Poisson's ratio of insect wings is unknown, we tested the effects of using a Poisson's ratio of 0.3 and found that the difference in model behavior was negligible. The uniform mass distribution of the models may cause dynamic bending in the distal regions to be overestimated, but qualitative differences in wing bending between the two models should not be affected. To determine the minimum number of elements necessary to capture the bending behavior of wings, we performed a sensitivity analysis with models composed of 200, 350, 865 and 2300 total elements, and found that 865 elements were sufficient to ensure asymptotic performance of the model.

We adjusted the material stiffness E of vein and membrane elements in the two models to produce different spatial patterns of flexural stiffness,but the same overall bending performance (so that tip and trailing edge displacement in response to a point load were the same as displacements measured in real wings). In the first model, we assigned all vein elements a single (homogeneous) Young's modulus of E=1.5×10⁸ N m^-2 (similar to values measured in locust hindwing; Smith et al., 2000), and all membrane elements a Young's modulus of E=2.1×10¹² N m^-2 (representing both increased material stiffness and the increased second moment of area of tubular veins; Fig. 2A). In the second model,we applied declining values of material stiffness to the model wing in 12 strips oriented diagonally (Fig. 2B); these strips are perpendicular to most of the wing veins,which decrease in diameter towards the wing edge and thus are likely to decrease in stiffness along this axis. We adjusted the values of material stiffness in these strips to approximate patterns of overall wing flexural stiffness measured in real Manduca wings (E in the model varies from 4.7×10⁷ N m^-2 to 4.5×10⁹ N m^-2 in membrane elements, and from 1.9×10¹¹ N m^-2 to 1.8×10¹³ N m^-2 in vein elements).

To determine the resulting pattern of flexural stiffness variation in the spanwise direction of the model wings, we fixed each wing at its base (with no displacement or rotation) and applied a point force at the tip. For chordwise measurements, we fixed the model at its leading edge (from the base to 2/3 span) and applied a point force at the trailing edge. We recorded the displacement at 22 nodes (junctions between elements) aligned between the point of attachment and the point of force application, and used this information to estimate spanwise and chordwise flexural stiffness distribution with the Matlab simplex minimization program, as in real wings.

Because the finite element program does not calculate aerodynamic forces acting on the wing, a damping factor had to be applied to stabilize the models. We compared the motions of the model wing with declining material stiffness to those of a real Manduca wing attached to a motor and rotated in the same way (Combes and Daniel,2003b; Daniel and Combes,2002). We found that a mass damping factor of 10 reduced high frequency vibrations and provided the closest match to motions of the real wing, and applied this damping factor to both model wings.

The simplex minimization analyses revealed that proposed EIdistributions based on either exponential or 2^nd-degree polynomial equations can provide relatively good predictions of measured wing displacement, while constant or linear distributions of EI cannot(see Combes, 2002). However,the proposed exponential EI distributions were more consistent within and between individuals, and effectively predicted displacement in the spanwise and chordwise directions of both Manduca sexta and Aeshna multicolor. Because these more consistent results allow us to make comparisons within and between individuals, we used an exponential equation to approximate flexural stiffness patterns in the species tested.

Our control measurements validated the procedure for estimating local flexural stiffness and revealed that the method is repeatable. Bending tests on a glass coverslip showed that a constant distribution of flexural stiffness provided the best approximation of measured displacement (as expected for a homogeneous beam). In addition, measurements of Manduca flexural stiffness in which each load was repeated several times provided consistent estimates of the EI distribution each time that a given load was applied (see Combes, 2002).

Morphological measurements showed that female Manduca sexta are significantly heavier and have larger wings than males (body mass, P=0.004; wing area, P=0.003; wing span, P=0.002,chord length, P=0.004; wing mass, P=0.012), but aspect ratio and wing loading were not significantly different in the individuals sampled(aspect ratio, P=0.946; wing loading, P=0.262). Although overall flexural stiffness scales with wing size across a broad range of species (see Combes and Daniel,2003a), spanwise

The wings of all Manduca sexta displayed dorsal/ventral asymmetry;the average flexural stiffness of their wings was higher in response to forces applied to the dorsal side of the wing, in both the spanwise and chordwise direction (spanwise, P=0.002; chordwise, P<0.001; Fig. 3A,B). In contrast, Aeshna multicolor displayed no dorsal/ventral difference in its hindwings or in the chordwise direction of its forewings (hindwing spanwise, P=0.327; hindwing chordwise, P=0.735; forewing chordwise, P=0.069; Fig. 3C,D). In the spanwise direction, the forewings were stiffer on the dorsal side than on the ventral side (P=0.018), but the difference between dorsal and ventral flexural stiffness was far smaller than that seen in Manduca.

The finite element model with diagonal strips of declining material stiffness accurately reproduced the sharp decline in flexural stiffness measured in male Manduca wings(Fig. 3A, red lines), while the model with homogeneous vein and membrane regions displayed a vastly different pattern of flexural stiffness (Fig. 3A, blue lines).

The effects of these differences in EI distribution were apparent in static tests on the model wings. Although an applied point force resulted in the same tip displacement, most bending in the exponential wing was confined to the outer third of the wing, while the homogeneous wing bent gradually along much of its span (Fig. 4A). When subjected to a pressure load, the homogeneous wing again bent along a large portion of its length, and its maximum displacement was nearly twice as large as the displacement of the exponential wing(Fig. 4B). In the exponential wing, displacement was localized to the tip and trailing edge of the wing, and bending was apparent in both the spanwise and chordwise directions.

The effects of an exponential decline in flexural stiffness were further illustrated by dynamic tests on the models. The homogeneous wing showed little chordwise bending and spanwise bending was most pronounced near the base,while the distal portion of the wing remained relatively rigid(Fig. 5, blue wings). In contrast, the exponential wing bent considerably in the chordwise direction,and bending in the spanwise direction was confined mainly to the outer portion of the wing (Fig. 5, red wings). For a movie of the model wings in motion, see http://faculty.washington.edu/danielt/movies.

Modeling passive wing flexibility requires some knowledge of the spatial distribution of stiffness throughout a wing. Although few studies have attempted to quantify spatial patterns of stiffness, a number of qualitative hypotheses have been put forward concerning the relationship between the structure of insect wings and their regional mechanical behavior. Most wings have relatively stiff, supporting zones near the wing base and leading edge,and relatively deformable areas near the edges of the wing(Wootton, 1981). Wing veins taper in diameter from base to tip(Wootton, 1992), and the presumed reduction in stiffness that this provides could serve many purposes,including reducing the inertial load at wing tips (thus reducing energy expenditure and stress at the wing base) and allowing wing tips to buckle and rebound from collisions without damage(Wootton, 1992). In addition,a wing that is stiff at the base and more flexible at the tip would be well suited to withstand the bending torques (which decline sharply from base to tip) generated by flapping (Ennos,1989).

Our measurements of spatial variation in flexural stiffness agree with these qualitative observations; flexural stiffness declines sharply from wing base to tip and from the leading to the trailing edge in the wings of both Manduca sexta and Aeshna multicolor. Displacement in the spanwise and chordwise directions of these wings can be predicted fairly well when the flexural stiffness distribution is approximated by an exponential equation. The observed similarities in flexural stiffness distribution in these two species (despite large differences in wing shape and venation pattern) suggest that a sharp decline in wing flexural stiffness towards the tip and trailing edge may be a common feature of many insect wings.

Although detailed inter-specific comparisons are not feasible with only two species, this study demonstrates several interesting intra-specific differences, particularly in Manduca sexta. Male and female Manduca sexta appear to differ in the spatial patterns of flexural stiffness in their wings; chordwise flexural stiffness declines far more sharply in male moths than in female moths(Fig. 3A,B), and as a result,the wings of female moths are significantly stiffer in the chordwise direction(dorsal

The wings of both male and female Manduca sexta display a large dorsal/ventral difference in average flexural stiffness, and in the exponents of flexural stiffness distribution. Average flexural stiffness is greater in response to forces applied on the dorsal side than on the ventral side (in both the spanwise and chordwise direction). In the spanwise direction, this difference is greatest near the base (Fig. 3A,B). Steppan(2000) found a similar dorsal/ventral asymmetry near the base of butterfly wings, and higher dorsal flexural stiffness has also been measured in the leading edge of locust hindwings (Wootton et al.,2000).

Although measurements of such bending asymmetry are rare, theoretical studies have suggested that insect wings might display dorsal/ventral stiffness asymmetry due to the camber inherent in most wings(Ennos, 1995, 1997; Wootton, 1993). We explored the effects of camber on bending asymmetry in the finite element model of a Manduca wing, adding 4% spanwise and 5% chordwise camber (the maximum values measured in real wings). Adding this degree of camber to the model wings had a relatively small effect on tip and trailing edge displacement (see Combes and Daniel, 2003a), and failed to produce any dorsal/ventral bending asymmetry (tip displacement was identical whether the force was applied to the convex or the concave side; see Combes, 2002).

Camber may in fact influence the maximum load each side of the wing can resist before buckling (the buckling load), and the finite element model does not simulate elastically stable buckling. However, our measurements on real wings do not display the strong non-linear dependence on loading that would be expected if elastically stable buckling were occurring during these static tests. In a more dynamic situation (such as flapping flight), where continuous shape changes can alter flexural stiffness, camber may in fact lead to dorsal/ventral differences in wing deformation. However, wing camber does not appear to explain the large differences measured in this study on wings subjected to static point loading.

The structural source of this bending asymmetry remains unclear. Although Manduca wings have none of the gross anatomical features often associated with bending asymmetry (e.g. a ventral flexion line), their veins may contain one-way hinges or other micro-structural features that facilitate asymmetric bending in other insect wings(Wootton, 1981; Wootton, 1992). In addition,the stress-stiffening effect explored by Kesel et al.(1998) in models of dipteran wings could potentially apply to Manduca. When unloaded, the membranes of Manduca wings lie rather loosely between the veins extending to the trailing edge. A point force applied to the dorsal (convex)side of the trailing edge appears to pull the trailing edge veins further apart, which would first remove any slack in the membrane, and then possibly increase the stiffness of the wing as force is applied to the already-taut membrane (the stress-stiffening effect). In contrast, when a point force is applied to the ventral side of the trailing edge, the veins are pushed together and the membranes between them become looser.

Although the finite element models do not reproduce this dorsal/ventral asymmetry, the strikingly different bending patterns of the two models demonstrate that flexural stiffness distribution is critical in determining how insect wings bend. In the model wing with exponentially declining flexural stiffness, both static and dynamic bending are limited to the tip and trailing edge, whereas bending in the homogeneous model occurs along the entire span.

These results suggest that the sharply declining flexural stiffness measured in real wings helps maintain rigidity near the wing base (despite larger bending moments), while localizing bending to the tip and trailing edge, which are regions of particular importance in controlling aerodynamic force production. The trailing edge is a critical control surface in airplane wings (affecting both the magnitude of lift and the lift-to-drag ratio; Anderson, 1991), and passive deformations in this region are likely to play a similarly important role in controlling flight forces in insects. The integration of passive flexibility and spatial patterns of flexural stiffness into models of insect flight will be an important step in determining the functional significance of wing structure and dynamic bending to insect flight performance.

If tip displacement is extremely large relative to beam length, the beam is bent down so far that its coordinates in the x-dimension while loaded are significantly different from its coordinates while unloaded, invalidating Equation A2. However, Equation A2 is valid for our measurements for the following reasons: (i) all measurements were standardized to a relative displacement (δ_T/L) of 0.05 in the spanwise direction and 0.08 in the chordwise direction, (ii) original measurements were rarely performed at relative displacements over 0.08 in the spanwise direction and 0.15 in the chordwise direction, and (iii) even in the most extreme cases(where relative displacement exceeded 0.15), the error in coordinates introduced by length changes in the x-dimension was less than 4%.

List of symbols
 
• a

  coefficient in exponential equation (A5)

 
• b

  intercept in linear equation (A4)

 
• c

  coefficient in exponential equation (A5)

 
• d constant

  (Equation A3)

 
• E

  material stiffness (Young's modulus)

 
• EI

  flexural stiffness

 
• \(EI\)

  average flexural stiffness

 
• f

  flapping frequency

 
• F

  force

 
• I

  second moment of area

 
• L

  beam length

 
• m

  slope in linear equation (A4)

 
• M

  moment

 
• p

  coefficient in polynomial equation (A6)

 
• q

  coefficient in polynomial equation (A6)

 
• r

  coefficient in polynomial equation (A6)

 
• t

  time

 
• x,y,z

  spatial coordinates

 
• x_o

  starting point of data collection

 
• x_max

  end point of data collection

 
• δ(x)

  displacement at position

 
• x δ_T

  displacement at tip or trailing edge

 
• δ_T/L

  relative displacement of tip or trailing edge

 
• θ

  rotation at wing base

 
• τ

  time constant

 
• ω

  circular frequency (2πf)

A. Trimble was critical in devising the method of measuring displacement by a laser ranging technique and in developing Matlab code to analyze loaded and unloaded images. M. Tu and E. Goldman generously provided useful comments and Matlab wisdom. J. Henry provided advice on collecting dragonflies and donated specimens to the study. J. Dierberger at MSC Software provided crucial troubleshooting of the FEM models, without which the dynamic simulations would not have been possible. This work was supported by NSF grant F094801 to T. Daniel, the John D. and Catherine T. MacArthur Foundation, an NSF graduate fellowship to S. Combes, and an ARCS fellowship to S. Combes.