SUMMARY

During flight, the wings of many insects undergo considerable shape changes in spanwise and chordwise directions. We determined the origin of spanwise wing deformation by combining measurements on segmental wing stiffness of the blowfly Calliphora vicina in the ventral and dorsal directions with numerical modelling of instantaneous aerodynamic and inertial forces within the stroke cycle using a two-dimensional unsteady blade elementary approach. We completed this approach by an experimental study on the wing's rotational axis during stroke reversal. The wing's local flexural stiffness ranges from 30 to 40 nN m2 near the root, whereas the distal wing parts are highly compliant (0.6 to 2.2 nN m2). Local bending moments during wing flapping peak near the wing root at the beginning of each half stroke due to both aerodynamic and inertial forces, producing a maximum wing tip deflection of up to 46 deg. Blowfly wings store up to 2.30 μJ elastic potential energy that converts into a mean wing deformation power of 27.3 μW. This value equates to approximately 5.9 and 2.3% of the inertial and aerodynamic power requirements for flight in this animal, respectively. Wing elasticity measurements suggest that approximately 20% or 0.46 μJ of elastic potential energy cannot be recovered within each half stroke. Local strain energy increases from tip to root, matching the distribution of the wing's elastic protein resilin, whereas local strain energy density varies little in the spanwise direction. This study demonstrates a source of mechanical energy loss in fly flight owing to spanwise wing bending at the stroke reversals, even in cases in which aerodynamic power exceeds inertial power. Despite lower stiffness estimates, our findings are widely consistent with previous stiffness measurements on insect wings but highlight the relationship between local flexural stiffness, wing deformation power and energy expenditure in flapping insect wings.

INTRODUCTION

For many years, the ability of insects to produce aerodynamic forces in excess of conventional steady-state aerodynamics has attracted considerable interest because it links aerial performance with the temporal distribution of lift and drag production (Dickinson et al., 1999). The majority of experimental studies using either robotic wings (e.g. Birch and Dickinson, 2003; Maybury and Lehmann, 2004; Sane and Dickinson, 2001) or computational approaches (e.g. Bos et al., 2008; Liu and Kawachi, 1998; Miller and Peskin, 2005; Ramamurti and Sandberg, 2002; Ramamurti and Sandberg, 2007; Sun and Lan, 2004; Wang et al., 2003) have thus focused on how insects employ unsteady aerodynamic mechanisms to enhance lift production during manoeuvring flight. Many previous studies on flapping insect wings have relied on the assumption that the wings of insects operate as rigid, non-deformable flat plates (Combes and Daniel, 2003a; Ho et al., 2003; Shyy, 2008; Vanella et al., 2009). However, similar to wing deformation in birds and bats, and unlike rigid plates, many insect wings undergo severe flexing and bending within a single stroke cycle, which may alter aerodynamic force production (Ennos, 1989; Rüppell, 1989).

Recently, Zhao et al. showed that as flexibility of robotic insect wings increases, its ability to generate forces monotonically decreases, whereas the lift-to-drag ratio remains approximately constant (Zhao et al., 2009). Their study further suggests that modulation of the trailing edge flexibility affects leading edge vorticity, thus providing a tool for manoeuvring control in insects. Numerical models on flapping flight in hover flies, by contrast, suggest that deformable wings produce approximately 10% more lift than rigid wings at approximately 5% reduced aerodynamic power requirements for flight (Du and Sun, 2010). Vanella et al. showed similar results for a chordwise-deforming two-dimensional wing blade (Vanella et al., 2009). According to their calculations, wing deformation increases the lift-to-drag ratio by up to 28% compared with a rigid blade. The recent work by Young and co-workers highlighted the impact of wing deformation on aerodynamic performance and power output in flying locusts (Young et al., 2009). Ennos was one of the first researchers to suggest that wing torsion at the stroke reversals of Diptera is primarily due to inertial forces occurring during the wings' acceleration and deceleration phases, and not to aerodynamic forces (Ennos, 1988; Ennos, 1989). However, robotic wing approaches and numerical modelling of these forces later demonstrated that lift and drag may peak during the stroke reversals and not during the translational part of the wing stroke (Dickinson et al., 1999). Combes and Daniel (Combes and Daniel, 2003a) showed that the prediction of the wing's dynamic shape change is thus essential for developing a comprehensive understanding of insect flight because instantaneous wing shape determines the direction and magnitude of fluid dynamic forces during wing flapping (Du and Sun, 2010; Shyy et al., 2008; Zhao et al., 2009).

Besides the various force components produced by flapping wings, wing deformation depends on the wing's primary structural properties: the spanwise stiffness and the chordwise flexibility. The leading edge of most insect wings is composed of a stiff structure with a three-dimensional relief, providing high rigidity in the spanwise direction (Combes and Daniel, 2003b). Thus, the flexural stiffness along the span, which is the overall bending stiffness of the wing, is one to two orders of magnitude greater than along the wing chord (Combes and Daniel, 2003b). Under aerodynamic loading, the wing's camber may change according to the wing's vein pattern, the veins' curvature and the joints between cross and longitudinal veins (Wootton, 1991). These three-dimensional skeletal networks of relatively stiff veins are interconnected through thin membranous and flexible cuticle cells (Hepburn and Chandler, 1976; Neville, 1960; Pfau, 1986). Because precise data about the local material properties of insect wings are rare, previous finite element approaches modelled the stiffness of the wing, i.e. the Young's modulus, equivalent to chitin, and simply used a stiffness constant for the entire wing surface (Kesel et al., 1998; Smith, 1996).

More recent studies have shown that wing stiffness can broadly vary between insect species and also throughout the wing surface, because elastic proteins such as resilin occupy key positions in the wing vein joints (Gorb, 1999; Haas et al., 2000a; Haas et al., 2000b). For example, the spatial variation in spanwise flexural stiffness covers approximately three orders of magnitude, ranging from approximately 10–6 N m2 in lacewings, Hemerobius, to 10–3 N m2 in the tobacco moth Manduca sexta. Steppan (Steppan, 2000) determined mean flexural stiffness over increasingly larger sections of dried Papilionoidea butterfly wings and found a flexural stiffness of 2.3×l0–8 to 1.49×10–6 N m2, and Wootton and co-workers (Wootton et al., 2000) measured deflections in three isolated sections of the locust rear wing. Combes and Daniel refined these measurements and determined the spatial pattern of flexural stiffness in 16 insect species from six orders by applying point forces on the wing in the spanwise and chordwise directions (Combes and Daniel, 2003b). In all tested species, spanwise flexural stiffness varied from 10–4 to 10–6 N m2 whereas chordwise flexural stiffness ranged between 10–5 and 10–7 N m2. The latter measurements also confirmed the dorso-ventral asymmetry in wing stiffness, as noted in previous studies (Ennos, 1988; Steppan, 2000; Wootton, 1993; Wootton et al., 2000).

In this study, we investigate the relative significance of aerodynamic and inertial forces during wing flapping for spanwise wing bending in the blowfly Calliphora vicina (Fig. 1). For this purpose, we measured wing mass and flexural stiffness of fly wings and numerically modelled the vertical and horizontal components of aerodynamic and inertial forces: lift and drag due to wing translation, wing rotation and Magnus force; and wing inertia due to the acceleration and deceleration of wing mass and added mass in each half stroke. The forces were derived by inserting a generic kinematic pattern for wing motion of flies and an experimental estimate for the wings' rotational axes into equations of a quasi-steady elementary blade approach published previously (Walker, 2002). We subsequently calculated and compared the displacement of discrete wing blades due to each force component and eventually evaluated to what degree potential energy stored during wing flexing may contribute to elastic energy recycling within the stroke cycle.

MATERIALS AND METHODS

Wing mass distribution

Inertial reaction forces during flapping motion depend on both wing acceleration and the distribution of wing mass. We estimated the spanwise mass distribution by cutting nine vein-dependent wing blade elements from freshly killed male Calliphora vicina Robineau-Desvoidy 1830 and using a UMX2 balance (Mettler Toledo, Inc., Columbus, OH, USA) with 0.1 μg precision. Thereafter, we converted these estimates into masses for blades with equal width (0.936 mm), yielding 93.0, 46.6, 23.8, 17.4, 18.3, 12.8, 11.0, 9.80 and 6.34 μg at distances from the wing root of 0.47, 1.40, 2.34, 3.27, 4.21, 5.21, 6.08, 7.02 and 7.95 mm, respectively (N=10 wings; Fig. 1K). The first moment of mass (normalized form) (Ellington, 1984a) was 43.5 (0.022), 65.3 (0.032), 55.8 (0.028), 57.0 (0.028), 77.1 (0.038), 65.6 (0.033), 66.7 (0.033), 68.8 (0.034) and 50.4 μg mm (0.025) for the nine blades. The centre of mass of each blade is most relevant for vertical inertia during wing rotation. We thus measured wing mass distribution chordwise, cutting the wings into seven stripes, parallel to the leading edge. Chordwise centre of mass is closely scattered around the rotational wing axis at a distance of approximately –48.1, –96.2, –144, –144, –96.2, –48.1, –96.2, –96.2 and 144 μm, respectively, where negative (positive) values indicate locations between the leading (trailing) wing edge and the rotational axis (Fig. 1K). The relative distance between the leading wing edge and the rotational axis and normalized to wing chord [i.e. xr in Walker (Walker, 2002)] is 0.25, 0.25, 0.23, 0.21, 0.21, 0.23, 0.24, 0.25 and 0.18 mm for the nine blades from wing base to tip. Because of the evaporation of hemolymph from the blade elements, we corrected the mass of all blades by a constant factor to match the sum of all blades to the mass of an intact wing. Total wing mass was approximately 239±67 μg (mean ± s.d., N=10 wings), and mean wing length and mean chord were approximately 8.42 and 3.01 mm, respectively. Mean wing area of each blade element was 1.81, 2.90, 3.16, 3.38, 3.38, 3.16, 2.81, 2.54 and 1.76 mm2, and the first moment of area (normalized form) (Ellington, 1984a) was 0.85 (0.004), 4.06 (0.019), 7.39 (0.035), 11.1 (0.053), 14.2 (0.068), 16.3 (0.078), 17.1 (0.082), 17.8 (0.085) and 14.0 mm3 (0.067), respectively.

Flexible stiffness measurements

To derive spanwise bending stiffness of the intact fly wing, we mounted the fresh wing with wax on a holder and applied a line load (metal wedge) to the wing via a glass spring at seven vein-dependent distances from the wing hinge (1.32, 2.51, 3.51, 4.70, 5.31, 6.11 and 7.30 mm; Fig. 1J). To avoid dry-out effects, the time between wing cutting and wax sealing was short and typically not longer than 1–2 min. Total experimental time was approximately 30 min. The strip load displaced the wing downward by up to 0.8 mm with linearly increasing magnitude, maintained the maximum deflection for 10 s and then returned to its initial position. The wedge was then advanced to a new spanwise position and another measurement was initiated. To avoid dynamic effects due to viscous changes of the cuticle, displacement velocity was set to approximately 130 μm s–1 (quasi-static conditions). While bending the wing, we recorded the displacement of both the wing and the glass spring using a fibre optic sensor (Fig. 1I). Reference displacements of the sensor conducted on a solid material (aluminium) were subsequently subtracted from the measurements and the spring constant was calculated from the slope between force and wing displacement. The spring constant of the measurement system was 150 N m–1 and thus at least three times the spring constant of the wing. Each experimental run produced up to 400 force and displacement measures and was conducted on both the ventral and dorsal wing surface, yielding a total of more than 6000 bending measures for four fly wings.

We validated our stiffness measurements by mounting a glass beam (19.7×0.99×0.14 mm) with wax to the experimental setup and conducted a series of measurements similar to the ones we performed on the fly wings. Data show that the wax did not hamper the compliance of the beam and did not cause any loss due to creep deformation. The spring constant of the glass beam was approximately 7 N m–1 and the Young modulus calculated from this measure was 68.7 GPa. This value is close to the known Young modulus of glass (72 GPa) and suggests that our measurements were within approximately 5% accuracy.

Fig. 1.

Wing deformation, methods and resilin distribution in Calliphora vicina. (A–H) Ventral wing flip in freely flying flies during take-off. The images were recorded using a Phantom V12 camera and show wing kinematics of two animals (A–D and E–H). Red lines and arrows highlight spanwise wing deformation in B and wing twisting in D. Image exposure time was 240 μs. Frames were recorded 50 ms (A,E) and 25 ms (B,F) prior ventral stroke reversal (C,G) and 25 ms after the reversal (D,H). (I) Experimental setup for measuring flexural stiffness in insect wings. The right wing from a freshly killed blowfly was attached to a holder and pressed downward with a metal wedge along seven chord lines (dotted red line) using a glass spring. Elastic force and displacement were measured via a fiber optic sensor. (J) Locations of flexural stiffness measurements. c, chord length. (K) Force calculations using a numerical approach were performaned on nine equally spaced wing blade elements. The centre of mass of each wing blade (red dot) is located near the rotational wing axis. w, width of wing blade; mw, wing mass of blade element; r, distance from wing root; R, wing length; RA, rotational axis during stroke reversal. (L,M) Wing areas containing the rubber-like protein resilin. Map of the resilin distribution (L, red) of right fly wings based on fluorescence microscopy of four flies. The area outlined in blue indicates the location of the resilin region shown in M.

Fig. 1.

Wing deformation, methods and resilin distribution in Calliphora vicina. (A–H) Ventral wing flip in freely flying flies during take-off. The images were recorded using a Phantom V12 camera and show wing kinematics of two animals (A–D and E–H). Red lines and arrows highlight spanwise wing deformation in B and wing twisting in D. Image exposure time was 240 μs. Frames were recorded 50 ms (A,E) and 25 ms (B,F) prior ventral stroke reversal (C,G) and 25 ms after the reversal (D,H). (I) Experimental setup for measuring flexural stiffness in insect wings. The right wing from a freshly killed blowfly was attached to a holder and pressed downward with a metal wedge along seven chord lines (dotted red line) using a glass spring. Elastic force and displacement were measured via a fiber optic sensor. (J) Locations of flexural stiffness measurements. c, chord length. (K) Force calculations using a numerical approach were performaned on nine equally spaced wing blade elements. The centre of mass of each wing blade (red dot) is located near the rotational wing axis. w, width of wing blade; mw, wing mass of blade element; r, distance from wing root; R, wing length; RA, rotational axis during stroke reversal. (L,M) Wing areas containing the rubber-like protein resilin. Map of the resilin distribution (L, red) of right fly wings based on fluorescence microscopy of four flies. The area outlined in blue indicates the location of the resilin region shown in M.

Wing elasticity and energy loss

Wing elasticity (elastic efficiency, ε) is the ratio between unloading (EU) and loading (EL) energy during wing bending:

This ratio is related to the resilience in rubber technology because blowfly wings contain the elastic protein resilin at various locations of the wing near the root (Andersen and Weis-Fogh, 1964) (Fig. 1L,M). In our stiffness measurements, energy was calculated from the product between elastic force and displacement (see Fig. 4). The wing's energy loss is the difference (1–ε) between maximum elasticity and the elasticity occurring during maximum wing deflection. The wing's elastic force decreases under static deflection and thus some of the elastic energy cannot be recovered under these experimental conditions (creep deformation; see Results).

formula
(1)

Wing kinematics

Stroke kinematic patterns in flies vary tremendously during manoeuvring flight because of steering muscle activity. In C. vicina, for example, both the activation of the sternobasalar muscle (M.b2) and phase advance of the first basalare (M.b1) is correlated with a rostral shift in downstroke kinematics, converting a figure-eight into an O-shaped wing tip trajectory pattern (Tu and Dickinson, 1996). Nachtigall reported similar kinematic alterations for the fly Phormia regina during tethered flight (Nachtigall, 1966). Mean stroke frequency in tethered C. vicina varies between 127 and 180 Hz (mean 158 Hz) (Nachtigall and Roth, 1983) and between 120 and 160 Hz (Balint and Dickinson, 2004), whereas slightly smaller values were found by Ennos for this species when flying freely (117 to 158 Hz) (Ennos, 1989). We confirmed these values by scoring stroke frequency within the first nine to 29 stroke cycles during unrestrained take-off behaviour, using a high-speed camera (Phantom V12, Vision Research, Wayne, NJ, USA). We found that frequency varied up to 30 Hz within a single fly and ranged from 131 to 201 Hz (means) between flies. Mean (±s.d.) frequency was 154±23.3 Hz (N=10 C. vicina). Stroke amplitude is also considerably variable in flies, ranging from 123 to 160 deg in C. vicina and from 136 to 180 deg in tethered and freely flying Drosophila (Ennos, 1989; Fry et al., 2003; Lehmann and Dickinson, 1998; Nachtigall and Roth, 1983; Zanker, 1990). It has further been shown that a delay in timing of wing rotation during the dorsal and ventral stroke reversal may attenuate force generation (Dickinson et al., 1999; Sane and Dickinson, 2002) and most insects such as beetles (Coccinella), dipterans (Calliphora, Tipula and Eristalis) and hymenopterans (Apis and Bombus) exhibit rather symmetrical wing rotation, in which 50% of the change in angle of attack occurs before and after each stroke reversal, respectively (Ellington, 1984b; Nachtigall, 1966).

We modelled wing kinematics in this study close to the patterns mentioned above and similar to the patterns used in other experimental and numerical studies on dipteran flight: 150 Hz stroke frequency, 135 deg flapping amplitude in a flat horizontal stroke plane, symmetrical wing rotation (equal rotational time during upstroke and downstroke), and an angle of attack at mid stroke of 40 and 20 deg during downstroke and upstroke, respectively (Dickinson et al., 1999; Lehmann and Pick, 2007) (Fig. 2). Using a rotational axis for wing rotation at 28.6% mean wing chord length (see Rotational wing axis), our kinematic pattern produces a mean vertical force opposite to gravity of 596 μN (two wings), supporting the fly's body weight within 2.8% accuracy (body mass, 59.1±10.2 mg, mean ±s.d., N=10 C. vicina). The same kinematic pattern produces 734 μN (two wings) with the rotational axis at the leading edge.

Rotational wing axis

Wing inertia and aerodynamic force production during the stroke reversals both depend on the location of the wing's rotational axis (see above). Previous studies often assumed that this axis runs from the wing base through the tip of the wing. Because of its significance for force production in this study, however, we estimated the rotational axis in tethered flying C. vicina more carefully. During stroke reversal, translational velocity of the wing is minimal but not zero, thus estimations of the rotational axis are inherently impaired by the wing's translational velocity. We coped with this limit for the analysis in the following way: we tethered blow flies to a wire using methods described previously (Lehmann and Dickinson, 1997), clipped the outer approximately 8% of the right wing and marked the wing at three positions using fluorescent dye (M1–M3, Fig. 3A). Wing motion was subsequently recorded at a capture rate of 6250 Hz with a Phantom V12 high speed video camera, orientated approximately normal to the longitudinal wing axis during the ventral stroke reversal (Fig. 3A,B). To achieve better temporal resolution, we numerically interpolated the positions of the three markers by a factor of approximately 14 and determined the stroke reversal by the minimum translational velocity of the wing tip (Fig. 3C). The wing's rotational axis was calculated by scoring the chordwise distances of two wing blades between the stroke reversal and time t(d) (downstroke, red, Fig. 3D) and between the stroke reversal and time t(u) (upstroke, black, Fig. 3D), with t(d) and t(u) 80 μs before and after the ventral reversal, respectively. The minimum distance indicates the rotational axis and is plotted for all tested animals (N=12 flies, 10 stroke cycles each; Fig. 3E). The difference in mean value between the downstroke (39.6% M1–M3) and upstroke estimates (50.9% M1–M3) reflects the contribution of wing translation to this analysis, causing overestimations during the upstroke and underestimations during the downstroke. The mean of all measures, converted into relative distance of mean wing chord length (, 3.01 mm), shows that the rotational axis is located 0.286 behind the leading edge.

Fig. 2.

Kinematic pattern used for numerical modelling. (A) Temporal changes of stroke angle (ϕ) for translational wing motion and geometrical angle of attack (α) with respect to the vertical. The kinematic pattern produces a flat wing tip trajectory in the horizontal with asymmetrical angle of attack during upstroke and downstroke and symmetrical wing rotation at the stroke reversals. (B,C) Absolute angular velocity and acceleration of the wing, respectively.

Fig. 2.

Kinematic pattern used for numerical modelling. (A) Temporal changes of stroke angle (ϕ) for translational wing motion and geometrical angle of attack (α) with respect to the vertical. The kinematic pattern produces a flat wing tip trajectory in the horizontal with asymmetrical angle of attack during upstroke and downstroke and symmetrical wing rotation at the stroke reversals. (B,C) Absolute angular velocity and acceleration of the wing, respectively.

Aerodynamic and inertial forces

We derived aerodynamic forces during wing flapping from a corrected version of a two-dimensional, semi-empirical unsteady blade element (USBE) model (Walker, 2002). The model covers all aerodynamic forces produced in a stroke cycle, including added mass reaction force and Magnus force. For fruit fly kinematics, the USBE model predicts 6–25% higher mean lift compared with computational fluid dynamic studies (Sun and Tang, 2002). Although presented throughout the manuscript, Magnus force is not considered for lift production and wing bending because it overestimates the rotational components during the stroke reversals compared with the forces measured in a three-dimensional robotic wing (Dickinson et al., 1999; Walker, 2002). The aerodynamic force coefficients required for USBE modelling were adopted from a data set measured in a three-dimensional electromechanical fly wing flapping at a Reynolds number of approximately 140 (Dickinson et al., 1999). Given the kinematic parameters mentioned above, mean Reynolds number for C. vicina is approximately 1268. We calculated instantaneous lift and drag including Magnus and added mass reaction forces due to fluid acceleration for each blade element and eventually integrated all forces spanwise to obtain total force during flapping motion (see the Appendix). Besides aerodynamic force, inertia is a major source for wing flexing and bending in insects (Combes and Daniel, 2003a; Ennos, 1988). These forces may simply be calculated from each wing blade mass element, translational and rotational wing motion (Maybury and Lehmann, 2004). Vertical inertia during in-plane flapping results from vertical displacement of wing mass during the rotational phases. Because the centre of mass for each blade is close to the wing's rotational axis, vertical inertia is small compared with horizontal inertia (Fig. 1K). The local bending moment of the fly wing eventually depends on the product between the moment arm, i.e. the distance of each blade element from the wing root, and the component of the vector sum of aerodynamic and inertial forces acting normal to the wing surface (see Results for illustration of angles). Thus, in a final step, we derived these components from a set of equations listed in the Appendix. Although the equations quantify wing deflection throughout the flapping cycle, we made no attempt to consider these changes for our calculations of aerodynamic forces (fluid structure interaction).

Fig. 3.

Experimental determination of the rotational wing axis in Calliphora vicina. (A) The right wing tip of intact flies was cut and marked (M1, M2, M3) using fluorescent dye. During wing flapping of tethered animals, marker positions were recorded before and after the ventral stroke reversal using high-speed video (images, right). (B) Typical example of wing movement at seven equally spaced times (0.16 ms sample period; t0 to t6) during the ventral stroke reversal. M1 (red) indicates the wing's leading edge and downstroke is from right to left. (C) Translational mean (±s.d.) wing tip velocity during stroke reversal. The stroke reversal occurs at minimum velocity. Times t(d) and t(u) indicate half sampling period before and after the ventral reversal, respectively. (D) Normalized chordwise distance (ΔS) between two wing positions at t(d) and stroke reversal (red), and between stroke reversal and t(u) (black). Minimum value indicates the wing's rotational axis. LE, leading wing edge; TE, trailing wing edge. (E) Location of rotational axis between the markers M1 and M3 for downstroke (red) and upstroke (black). Axis of wing rotation is approximately 28.6% mean wing chord length (3.01 mm, shown in A) behind the wing's leading edge. Data are means ± s.d., N=120 stroke cycles from 12 flies.

Fig. 3.

Experimental determination of the rotational wing axis in Calliphora vicina. (A) The right wing tip of intact flies was cut and marked (M1, M2, M3) using fluorescent dye. During wing flapping of tethered animals, marker positions were recorded before and after the ventral stroke reversal using high-speed video (images, right). (B) Typical example of wing movement at seven equally spaced times (0.16 ms sample period; t0 to t6) during the ventral stroke reversal. M1 (red) indicates the wing's leading edge and downstroke is from right to left. (C) Translational mean (±s.d.) wing tip velocity during stroke reversal. The stroke reversal occurs at minimum velocity. Times t(d) and t(u) indicate half sampling period before and after the ventral reversal, respectively. (D) Normalized chordwise distance (ΔS) between two wing positions at t(d) and stroke reversal (red), and between stroke reversal and t(u) (black). Minimum value indicates the wing's rotational axis. LE, leading wing edge; TE, trailing wing edge. (E) Location of rotational axis between the markers M1 and M3 for downstroke (red) and upstroke (black). Axis of wing rotation is approximately 28.6% mean wing chord length (3.01 mm, shown in A) behind the wing's leading edge. Data are means ± s.d., N=120 stroke cycles from 12 flies.

Elastic potential energy and strain energy density

Elastic potential energy is the energy stored as a result of deflection of the elastic wing and was approximated by the work done to deflect a cantilever beam. We calculated the local elastic potential energy (EPE) from local deflection (δ) and local flexural stiffness (EI) of each wing blade at a given distance (r) from the wing root using the following equation:

where w is blade width, and subsequently integrated the local value to determine total elastic energy within the stroke cycle. The local strain energy density (SED) for each wing blade equals the work stored in the wing during bending divided by the volume of each wing blade. We calculated this measure from the ratio between local bending moment (M), EI, c and wing thickness (h), given by:

formula
(2)

Wing thickness was derived from total wing mass and area, assuming a density for chitin of 1200 kg m–3 (Ellington, 1984a). For the nine blades from root to tip we obtained thicknesses of 42.8, 13.4, 6.28, 4.30, 4.52, 3.36, 3.25, 3.22 and 3.00 μm, respectively (Fig. 1K). In contrast to EPE, we averaged SED spanwise to obtain mean SED during wing flapping.

formula
(3)

Fluorescence microscopy of resilin

Resilin is a highly elastic protein occurring in various parts of the flight apparatus of the fly. We determined the distribution of resilin in wings freshly cut from C. vicina, mounting them on coverslips in a water-soluble medium (Moviol, Hoechst, Frankfurt, Germany). Resilin has an autofluorescence at a very narrow band of wavelengths of approximately 400 nm. Thus, resilin in biologically native structures can be determined without immune labelling (Andersen and Weis-Fogh, 1964). To enhance resilin visibility, we observed the fly wings at three wavelengths of a fluorescent microscope (Zeiss Axioskop, Zeiss, Oberkochen, Germany): green (excitation 512–546 nm, emission 600–640 nm), red (excitation 710–775 nm, emission 810–890 nm) and ultraviolet (excitation 340–380 nm, emission 425 nm). We recorded images in each of the three bands of wavelengths using a video camera (Sony 3CCD camera DXC-950P) and superimposed them as semi-transparent pictures using Adobe Photoshop (Fig. 1L,M).

RESULTS

Elasticity, spring constant and local flexural stiffness

The elastic force of blowfly wings linearly increased with increasing load applied chordwise until the displacement was constant (t0, Fig. 4A). Under constant load, the elastic force typically decreased with increasing measurement time due to creep deformation, which produces a loss in energy due to force and displacement (exponential decay; upward bending, 8.82 s; downward bending, 6.02 s; see legend of Fig. 4B,D,E). Creep deformation was close to saturation at the beginning of the unloading procedure (offset of fit curve; upstroke, 0.828 mN; downstroke, 1.52 mN; Fig. 4D,E). Elasticity for seven wing blades is plotted in Fig. 4F and shows a maximum of 20–23% loss of elastic energy due to creep deformation under constant bending stress. Mean elasticity was not significantly different between upward and downward bending (t-test, upward, 0.77±0.14, N=24; downward, 0.80±0.10, N=15; P=0.45). The same holds for the blade-wise comparison between both bending directions (ANOVA, blade at 1.32 mm, P=0.31; 2.51 mm, P=0.45; 3.51 mm, P=0.83; 4.70 mm, P=0.88; 5.31 mm, P=0.53; 6.16 mm, P=0.27). Moreover, the elasticity of the wing blades does not significantly increase with increasing distance from the wing root (linear regression, y=0.82–7.91x, R2=0.01, P=0.51, N=39; Fig. 4F).

Fig. 4.

Stiffness and elasticity measurements of Calliphora vicina wings. (A) Change in elastic force (black) and combined displacement of sensor and wing (red) during wing loading experiments. Dotted lines indicate range of constant displacement. (B) Elastic force plotted against combined displacement of sensor and wing for two distances from the wing root. (C) Loading and unloading area during wing displacement as shown for the example in B (black). These data were used to calculate elasticity (shown in F). (D,E) Elastic force during constant wing bending in the dorsal and ventral directions, respectively (10 s). Red line shows exponential fit (D, y=828+142ex/8.82, χ2/DoF=542, R2=0.74, N=3424; E, y=1518+187ex/6.07, χ2/DoF=296, R2=0.86, N=2747; data of each blade are normalized to the same mean). (F) Elasticity indicating normalized loss of elastic force during bending in the dorsal (open circles) and ventral directions (filled circles). (G) Spring constants averaged from root to the measured point derived from wing bending in the dorsal (white bars) and ventral directions (grey bars), with red circles indicating the means of both. (H) Combined flexural stiffness and polynomial fit (red, y=–49.5+71.1x–7.28x2, R2=0.38) calculated from data shown in G. (I) Local flexural stiffness of fly wings and the logistic fit function {red, y=33.0[1+e1.70(x–4.10)]–1, χ2/DoF=9.05, R2=0.96}. Data are means ± s.d.

Fig. 4.

Stiffness and elasticity measurements of Calliphora vicina wings. (A) Change in elastic force (black) and combined displacement of sensor and wing (red) during wing loading experiments. Dotted lines indicate range of constant displacement. (B) Elastic force plotted against combined displacement of sensor and wing for two distances from the wing root. (C) Loading and unloading area during wing displacement as shown for the example in B (black). These data were used to calculate elasticity (shown in F). (D,E) Elastic force during constant wing bending in the dorsal and ventral directions, respectively (10 s). Red line shows exponential fit (D, y=828+142ex/8.82, χ2/DoF=542, R2=0.74, N=3424; E, y=1518+187ex/6.07, χ2/DoF=296, R2=0.86, N=2747; data of each blade are normalized to the same mean). (F) Elasticity indicating normalized loss of elastic force during bending in the dorsal (open circles) and ventral directions (filled circles). (G) Spring constants averaged from root to the measured point derived from wing bending in the dorsal (white bars) and ventral directions (grey bars), with red circles indicating the means of both. (H) Combined flexural stiffness and polynomial fit (red, y=–49.5+71.1x–7.28x2, R2=0.38) calculated from data shown in G. (I) Local flexural stiffness of fly wings and the logistic fit function {red, y=33.0[1+e1.70(x–4.10)]–1, χ2/DoF=9.05, R2=0.96}. Data are means ± s.d.

The combined spanwise spring constant decreases with increasing distance from the wing root (Fig. 4G) whereas the combined flexural stiffness of 124×10–9 N m2 saturates at approximately 4.8 mm or 43% wing length. For comparison, this value is approximately eight times lower than a value previously measured for Calliphora sp. at 70% wing length (Combes and Daniel, 2003b). None of the spring constants between upward and downward bending were significantly different in each blade element (t-test, P>0.05 for all blade elements). Because of thicker veins, the wing root yields an elevated mean local flexural stiffness of approximately 33.1×10–9 N m2, whereas the more distal parts of the wing are relatively compliant, exhibiting a local flexural stiffness ranging from 0.64×10–9 to 2.29×10–9 N m2 (Fig. 4I).

Aerodynamic and inertial forces

The USBE model allowed us to simulate the instantaneous force components that potentially affect wing bending in flapping C. vicina wings: the components of aerodynamic forces including added mass reactions forces (Fig. 5A), inertial forces due to acceleration and deceleration (Fig. 5B), and Magnus force, normal to the wing surface (Fig. 5C). Fig. 5D–G shows how these forces vary throughout the stroke cycle. Magnus force points upwards, pulling on the wing in the direction of lift, whereas inertia mainly acts in the horizontal (Fig. 5D,G). Horizontal wing motion due to translation contributes approximately 33 times more to maximum horizontal inertia than horizontal wing motion due to wing rotation (peak forces; Fi,h*=1.92 mN, Fi,h**=0.06 mN, respectively; Fig. 5F) because wing mass is closely scattered around the wing's longitudinal axis, and off-axis masses are thus low (Fig. 1K). Vertical inertia due to wing rotation is always negative and peaks at small values of approximately –0.07 mN (red, Fig. 5F). This is of interest, as the conventional treatment of inertial power ignores the contribution of horizontal and vertical inertia due to wing rotation (Ellington, 1984c). Averaged throughout the cycle, total absolute aerodynamic force (Fa,n) and inertia (Fi,n) normal to the wing surface amount to 474 and 235 μN, respectively, whereas added mass reaction force (166 μN) as part of Fa,n amounts to 35.0 and 70.6% of Fa,n and Fi,n, respectively. Similar values were obtained for the component of mean Magnus force acting normal to the wing surface (Fm,n=56.8 μN; 12.0% of mean Fa,n; 24.2% of mean Fi,n).

Fig. 5.

Forces acting on the center of wing mass mw for each blade element during flapping motion. (A) Total aerodynamic force Fa is the vector sum of lift (L), drag (D) and added mass reaction force (Facc). (B) Total inertial force Fi depends on horizontal (Fi,h) and vertical inertia (Fi,v) due to wing translation (*) and rotation (**), respectively. (C) Magnus force Fm points vertically upward. The force components normal to the wings surface Fa,n (blue), Fi,n (green) and Fm,n (red) cause wing bending throughout the stroke cycle. For A–C, black lines indicate a wing blade element; black circles indicate the wing's leading edge. αg, geometrical angle of attack with respect to horizontal; i, three-quarter chord location of attached flow vector (Walker, 2002); RA, rotational axis of wing. (D) Alterations in forces normal to the wing surface during downstroke (top) and upstroke (bottom). d, dorsal; v, ventral. (E) Vertical (blue, Fa,v) and horizontal (red, Fa,h) aerodynamic forces, both including Facc, derived from lift, drag and angle of local stream [see eqns 9 and 10 in Walker (Walker, 2002)]. (F) Wing-mass-dependent vertical and horizontal inertia and added mass reaction force due to inertial forces by the fluid Facc. (G) Aerodynamic, inertial and Magnus force components normal to the wing surface. Magnus force is not considered for wing bending analysis (see Materials and methods). Forces shown in D–G are total values, combined spanwise among all wing blade elements.

Fig. 5.

Forces acting on the center of wing mass mw for each blade element during flapping motion. (A) Total aerodynamic force Fa is the vector sum of lift (L), drag (D) and added mass reaction force (Facc). (B) Total inertial force Fi depends on horizontal (Fi,h) and vertical inertia (Fi,v) due to wing translation (*) and rotation (**), respectively. (C) Magnus force Fm points vertically upward. The force components normal to the wings surface Fa,n (blue), Fi,n (green) and Fm,n (red) cause wing bending throughout the stroke cycle. For A–C, black lines indicate a wing blade element; black circles indicate the wing's leading edge. αg, geometrical angle of attack with respect to horizontal; i, three-quarter chord location of attached flow vector (Walker, 2002); RA, rotational axis of wing. (D) Alterations in forces normal to the wing surface during downstroke (top) and upstroke (bottom). d, dorsal; v, ventral. (E) Vertical (blue, Fa,v) and horizontal (red, Fa,h) aerodynamic forces, both including Facc, derived from lift, drag and angle of local stream [see eqns 9 and 10 in Walker (Walker, 2002)]. (F) Wing-mass-dependent vertical and horizontal inertia and added mass reaction force due to inertial forces by the fluid Facc. (G) Aerodynamic, inertial and Magnus force components normal to the wing surface. Magnus force is not considered for wing bending analysis (see Materials and methods). Forces shown in D–G are total values, combined spanwise among all wing blade elements.

At the stroke reversals, inertial peaks due to horizontal translational deceleration and acceleration clearly dominate total bending force normal to wing chord, outscoring aerodynamic force approximately 2.24-fold (maximum Fi,n=1.92 mN; maximum Fa,n=0.86 mN). However, the normal components of wing inertia and lift/drag run counter at the beginning of the stroke reversal and thus partly cancel each other out (Fig. 5G). Normal Magnus force shows an elaborated time course but is too small to cause significant bending (Fig. 5G; see Materials and methods). Superficially, the data suggest that maximum bending force in the blowfly is primarily caused by inertial effects at the stroke reversals and not by aerodynamic force production. Changes in the shape of the translational velocity profile during back-and-forth flapping motion should thus have a predominant effect on the magnitude of wing bending in flies. However, spanwise wing deformation does not depend on local forces but on local bending moments and thus on how the forces are distributed in the spanwise direction during wing flapping.

Local bending moments

Local bending moments depend on the spanwise distribution of normal-orientated aerodynamic, inertial and Magnus forces that vary distinctly throughout wing flapping (Fig. 6). Because local forces and local bending moments are rather similar during upstroke and downstroke, Fig. 6B–H only shows data for the downstroke. As already mentioned, at the beginning and the end of each half stroke, during the wing acceleration and deceleration phases, respectively, inertial forces peak and cause a homogenously distributed spanwise load on the wing (time bin 0–0.1, Fi,n=87.9±16.2 μN; time bin 0.4–0.5, Fi,n=–49.6±7.70 μN; means ± s.d., N=9 wing blades, downstroke; Fig. 6B,F). This effect results from the combined changes in translational and rotational velocities, although the contribution of rotational velocities is relatively small (see above, Fig. 1K and Fig. 5F). In contrast to the stroke reversals, inertial forces are low during wing translational phases and thus do not contribute to wing bending moments (Fig. 6C–E). Aerodynamic force, by contrast, increases from root to tip due to the spanwise increase in translational velocity of the local blades (time bin 0–0.1, Fa,n= 49.4 ±27.3 μN; time bin 0.4–0.5, Fa,n=114±56.7 μN; means ± s.d., N=9 wing blades, downstroke), producing only minor fluctuations in local bending moments over time (Fig. 6B,E). Aerodynamic-induced bending, however, peaks at the end of each half stroke when translational and rotational wing velocity is maximum (Fig. 6F). As already shown in Fig. 5, Magnus force, is relatively small and contributes, if considered, moderately to local bending (time bin 0–0.1, Fm,n=–8.06±4.45 μN; time bin 0.4–0.5, Fm,n=8.65±4.78 μN; means ± s.d., N=9 wing blades, downstroke; Fig. 6B,F).

Fig. 6.

Spanwise force on wing blades and local bending moments. (A) According to cantilever beam theory, local bending moments are calculated from the product between moment arms and forces acting distal to the local blade. (B–F) Aerodynamic (blue, Fa,n) and inertial force (red, Fi,n, right scale) orientated normal to the wing surface. Local bending moments (line, Mi, Ma) are plotted in the corresponding color with moments caused by Magnus force (Mm) in grey. Data are averaged within binned normalized time of the stroke cycle lasting from 0 to 1 (B, 0–0.1; C, 0.1–0.2; D, 0.2–0.3; E, 0.3–0.4; F, 0.4–0.5). (G) Spanwise total forces (time bins 1–5 as in B–F). Data in B–G are shown for the downstroke. (H,I) Spanwise total bending moments during downstroke in H (time bins as in G) and upstroke in I (time bins 1–5 are 0.5–0.6, 0.6–0.7, 0.7–0.8, 0.8–0.9 and 0.9–1.0 normalized time, respectively).

Fig. 6.

Spanwise force on wing blades and local bending moments. (A) According to cantilever beam theory, local bending moments are calculated from the product between moment arms and forces acting distal to the local blade. (B–F) Aerodynamic (blue, Fa,n) and inertial force (red, Fi,n, right scale) orientated normal to the wing surface. Local bending moments (line, Mi, Ma) are plotted in the corresponding color with moments caused by Magnus force (Mm) in grey. Data are averaged within binned normalized time of the stroke cycle lasting from 0 to 1 (B, 0–0.1; C, 0.1–0.2; D, 0.2–0.3; E, 0.3–0.4; F, 0.4–0.5). (G) Spanwise total forces (time bins 1–5 as in B–F). Data in B–G are shown for the downstroke. (H,I) Spanwise total bending moments during downstroke in H (time bins as in G) and upstroke in I (time bins 1–5 are 0.5–0.6, 0.6–0.7, 0.7–0.8, 0.8–0.9 and 0.9–1.0 normalized time, respectively).

During the wing acceleration phase, aerodynamic/added mass reaction force and inertia produce local bending moments of up to approximately 5.81 μN m (downstroke) and –4.71 μN m (upstroke) at the wing root, whereas during wing deceleration bending is significantly smaller (downstroke, 3.41 μN m; upstroke, –2.84 μN m; Fig. 6G–I). Wing bending thus occurs predominantly at the beginning of each half stroke, because at the end of each half stroke aerodynamic and inertial bending moments act in opposite directions, pulling and pushing on the dorsal (ventral) wing surface during the downstroke (upstroke), respectively. At this time, aerodynamic force is apparently slightly more effective for spanwise wing bending because of the pairing between elevated force and the long moment arm (Fig. 6F). The force balance also suggests that inertial power is put straight into aerodynamic power instead of wing deflection power and so elasticity is not required to limit inertial power requirements at the end of each half stroke (Ellington, 1984c). The sign changes of spanwise total force further produce negative and positive bending moments on proximal and distal wing blades at the stroke reversals as shown in Fig. 6G (time bin 0.4–0.5; negative bending, 0–1.7 mm wing length). If the chordwise positions of total force vectors acting normal to the wing surface differ slightly between blades, opposing bending moments promote wing twisting. Surprisingly, we observed wing twisting near the wing root 25 ms after the ventral stroke reversal in freely flying C. vicina (arrow, Fig. 1D).

Wing tip deflection, EPE and SED

Local wing deflection depends on the ratio between local bending moments and local flexural stiffness (see Eqn A10 in the Appendix). Consequently, local deflection is elevated near the wing root because of high bending moments and despite high local stiffness, and also elevated near the wing tip due to small local stiffness and despite small bending moments (Fig. 7E). During the wing acceleration phases, local bending is maximum and amounts to 61.7 μm at the root and 84.7 μm for a blade near the wing tip, whereas the wing's middle parts are less prone to deflection (Fig. 7E). The temporal change in wing tip deflection for two stroke cycles and as calculated from Eqn A11 is shown in Fig. 7B. As already suggested, this measure indicates that spanwise deformation peaks at the beginning of the half stroke (dorsal reversal, 46.2 deg; ventral reversal, 43.4 deg); the wing is less deformed throughout the upstroke and downstroke (Fig. 7A, middle and right). We found comparable wing deformations (approximately 27 deg) at the beginning of the stroke reversal in freely flying C. vicina using high-speed video recording (Fig. 7A, left).

Fig. 7.

Wing tip deflection (δ), elastic potential energy (EPE) and strain energy density (SED) during wing flapping. (A) Left: wing bending at the beginning of the upstroke in a freely flying Calliphora vicina. Body and wing shape were redrawn from the high-speed image shown in Fig. 1. Middle and right: spanwise wing bending for 20 equally spaced times (0.05 stroke cycle) during upstroke and downstroke based on wing kinematics, stiffness measurements and numerical force modelling. (B) Wing tip deflection. (C) Total EPE stored by the entire wing during wing flapping and (D) the wing's mean SED. Red line shows smoothed data (B-spline on 40 equally spaced data bins). (E) Deflection of the distal end of the blade, (F) local EPE and (G) local SED during downstroke. Data are shown for five normalized time bins of the stroke cycle (red, 0–0.1; black, 1, 0.1–0.2; 2, 0.2–0.3; 3, 0.3–0.4; and blue, 0.4–0.5 normalized time). See legend of Fig. 6 for more information.

Fig. 7.

Wing tip deflection (δ), elastic potential energy (EPE) and strain energy density (SED) during wing flapping. (A) Left: wing bending at the beginning of the upstroke in a freely flying Calliphora vicina. Body and wing shape were redrawn from the high-speed image shown in Fig. 1. Middle and right: spanwise wing bending for 20 equally spaced times (0.05 stroke cycle) during upstroke and downstroke based on wing kinematics, stiffness measurements and numerical force modelling. (B) Wing tip deflection. (C) Total EPE stored by the entire wing during wing flapping and (D) the wing's mean SED. Red line shows smoothed data (B-spline on 40 equally spaced data bins). (E) Deflection of the distal end of the blade, (F) local EPE and (G) local SED during downstroke. Data are shown for five normalized time bins of the stroke cycle (red, 0–0.1; black, 1, 0.1–0.2; 2, 0.2–0.3; 3, 0.3–0.4; and blue, 0.4–0.5 normalized time). See legend of Fig. 6 for more information.

Because of high stiffness near the wing root and elevated bending moments, most of the EPE during spanwise bending is stored in the proximal parts of the wing (Fig. 7F). Maximum EPE of a local blade amounts to 0.19 μJ at the beginning of the half stroke (red, Fig. 7F), which is almost completely released during the wing's translation phase. By integrating EPE spanwise, we estimated total EPE as shown in Fig. 7C and calculated mean EPE for a full half stroke (downstroke, mean EPE=0.27 μJ; upstroke, mean EPE=0.18 μJ). Local SED is a measure of how much strain energy is stored during spanwise deformation per unit wing volume. In contrast to EPE, this measure does not significantly change in the spanwise direction, despite the massive change in local wing volume (see Materials and methods). Bending stress is thus concentrated on the proximal wing part (0–0.1 stroke cycle, 6.92±1.59×103 J m–3; 0.4–0.5 stroke cycle, 2.81±0.55×103 J m–3; SED mean ± s.d., N=8 blades, downstroke; Fig. 7G). In this respect, the distribution of resilin at the proximal end of the wing is of particular interest. Although elasticity does not necessarily increase the toughness of the wing, resilin patches might reduce the risk of breaking near the wing hinge because of a reduction in peak stress on the rigid parts of the wing during flapping (Fig. 1L,M, Fig. 7D).

DISCUSSION

The aim of this study was to evaluate the relative contribution of aerodynamic and inertial forces to wing bending in the blowfly C. vicina and to estimate how much elastic potential energy is stored in a fly wing during flapping motion. Because we restricted our modelling to cantilever analysis of spanwise deformations, our data cannot explain complicated chordwise deflection or bending due to wing vein hinges (Wootton, 1992). The combination of USBE modelling, a numerical framework for wing inertia, and stiffness measurements on fly wings revealed that fly wings might bend spanwise more than 40 deg at the beginning of each half stroke (Fig. 7). We observed smaller but still comparable deformation in C. vicina taking off freely from a start platform (Figs 1 and 7). Bending moments at the stroke reversals are mainly driven by two forces: horizontal inertia due to the translational velocity profile of the wing and aerodynamic forces due to wing rotation. Both forces produce an increasing local bending moment with decreasing distance from the wing root, whereas the high local flexural stiffness at the wing root avoids tremendous deformations of the proximal blades (Fig. 6). In turn, despite small bending moments, wing deflection of distal wing blades is significant because of the increase in structural compliance from root to tip (Fig. 7).

Wing stiffness

The exponential decay in wing stiffness as shown in Fig. 4G is consistent with previous measurements on insect wings, despite the different measurement techniques [line vs point measurements in Combes and Daniel (Combes and Daniel, 2003b; Combes and Daniel, 2003c)] and bending methods [segmental force application vs single force application at 70% wing length in Combes and Daniel (Combes and Daniel, 2003b; Combes and Daniel, 2003c)]. However, the combined flexural stiffness at 70% wing length of Calliphora sp. in Combes and Daniel (Combes and Daniel, 2003b; Combes and Daniel, 2003c) is approximately eight times higher (approximately 10–6 N m2) than our measurements (1.2×10–7 N m2; Fig. 4H). The origin of this difference is unclear but might result from dry-out effects in the Combes and Daniel studies. For comparison, we calculated spanwise bending using the higher stiffness value by Combes and Daniel. According to our aerodynamic modelling, the higher stiffness estimate would limit wing-tip bending to a few degrees, making the wing virtually rigid during wing flapping and running counter to the spanwise deformations derived from freely flying blowflies (Figs 1 and 7).

The high flexural stiffness near the wing root is consistent with studies on vein topology (Wootton, 1991; Vanella et al., 2009). In contrast to the low compliance of the wing root, Wootton suggested that the high compliance of distal wing parts might serve as a mechanism for wing protection when the wing hits solid objects during flapping motion (Wootton, 1992). In this case, high elasticity at the tip would allow the tip to buckle and rebound without further damage (Combes and Daniel, 2003c). Wing-tip damage attenuates aerodynamic performance, yielding two effects: an increase in wing stroke frequency due to a change in load of the thoracic oscillator and, potentially, a reduction in the ability to generate a stable leading edge vortex due to changes in tip vorticity [fruit fly (Bender and Dickinson, 2006); bumblebee (Hedenström et al., 2001)].

Although the asymmetry for upward and downward bending as shown in Fig. 4D,E apparently supports previous work on wings of butterflies, the hawkmoth Manduca sexta, and locust hindwings (Combes and Daniel, 2003b; Steppan, 2000; Wootton et al., 2000), its structural origin remains unclear, as pointed out by Combes and Daniel. A possible explanation of this phenomenon might be the stress-stiffening effect of dipteran wings (Kesel et al., 1998). According to previous suggestions, stress stiffening may result from the intrinsic motion of veins and membranes under load (Combes and Daniel, 2003b). A force pressing on the dorsal wing side might pull the trailing edge veins further apart, which removes slack in the membrane, whereas forces on the ventral side of the wing push the veins together and reduce membrane tension.

Energy storage and elasticity

It has previously been shown that elastic energy recycling may limit the energetic expenditure during flapping flight of insects (Casey and Ellington, 1989; Dickinson and Lighton, 1995; Ellington, 1984c). In the fruit fly Drosophila, for example, a minimum of 10% elastic storage is required to minimize flight costs during hovering flight conditions (Dickinson and Lighton, 1995). However, this calculation ignored wing drag due to three-dimensional unsteady aerodynamic phenomena, thus the need for elastic storage of kinetic energy during hovering is presumably less than originally suggested (Lehmann, 2001). Nevertheless, the storage of kinetic energy due to elastic wing bending might contribute to the overall energy balance during wing flapping. We thus found it surprising that the wing does not gain much energy during its deceleration phases, because the moments produced by lift/drag and added mass reaction force partly cancel out wing mass-induced moments (Fig. 6). Total EPE at the end of the upstroke and downstroke (0.9–1.0 and 0.4–0.5 relative time, respectively) is thus relatively small, amounting to only 0.43 and 0.64 μJ, respectively (Fig. 7C). Consequently, the elastic structures of the wing do not recycle much kinetic energy gained from a preceding half stroke and thus contribute only slightly to the recycling of kinetic energy at the stroke reversals.

A higher amount of EPE is stored, by contrast, at the beginning of each half stroke, up to approximately 1.89 μJ (upstroke) and 2.30 μJ (downstroke), and subsequently released throughout the wing translation phase. These measurements equate to approximately 2.5 (upstroke) and 3.1 (downstroke) times the mean kinetic energy required for wing acceleration in the first half of each half stroke (upstroke, 0.76 μJ; downstroke, 0.74 μJ) (Ellington, 1984c). Because wing bending delays the sudden acceleration of local wing mass, spanwise wing deformation might lead to a temporal reduction in inertial costs at times of maximum energetic expenditures (i.e. wing acceleration phases) within the stroke cycle. Benefits of elasticity for force production and muscle dynamics have previously been described in locomotion (e.g. Biewener and Roberts, 2000; Lichtwark and Barclay, 2010; Roberts et al., 1997). If we assume that elastic deformation of the wing helps to distribute instantaneous flight power requirements more homogenously within the stroke cycle, this mechanism could lower the peak requirements of flight muscle force production. Further considering that flight muscle mass is largely determined by the need to produce peak mechanical power, elastic wing deformation might help to reduce flight muscle size, in turn increasing the payload capacity of the flying insect.

Our stiffness measurements, however, suggest that not all potential energy stored during wing bending is necessarily released. By estimation of elasticity during constant wing loading, we found a loss of energy ranging from 20 to 23%, or 77 to 80% elastic efficiency (Fig. 4F). Previous estimates of elastic efficiency of the locust thorax at normal work conditions indicate similar values: intact thorax, 86%; empty box, 88%; elastic hinges, 97%; and non-active muscle, 80% (Andersen and Weis-Fogh, 1964). A more recent study reports 95% elastic efficiency for resilin (Lv et al., 2010). Because elastic deformation occurring under static loading certainly differs from that occurring during high-frequency loading and unloading cycling in naturally flapping fly wings, it is likely that the obtained value presents a lower rather than an upper estimate for elastic loss. A study on the loss factor (1–elastic efficiency) of an elastic tendon of the locust, for example, suggests that inelasticity increases with increasing frequency of alternating strain from approximately 1% at 25 Hz to approximately 5% at 100 Hz (Jensen and Weis-Fogh, 1962). Eventually, mean wing deformation power is the product between elastic efficiency (0.8), mean EPE (228 nJ) and stroke frequency (150 Hz), and amounts, in C. vicina, to 27.3×10–3 mW, averaged over the entire stroke cycle. This value is relatively small and equates to approximately 5.9% of the inertial power requirements for wing acceleration (upstroke, 0.46 mW; downstroke, 0.47 mW) and approximately 2.3% of aerodynamic power requirements due to horizontal drag (upstroke, 0.81 mW; downstroke 1.58 mW; see Eqns A13A14 in the Appendix). These data suggest that the loss in wing deformation power due to wing elasticity is relatively small (6.83×10–3 mW).

CONCLUSIONS

In conclusion, investigations of flapping wing deformation in insects fuel research on the interference between flight forces and wing structures, termed fluid–structure interactions (Shyy et al., 2008). Even small deformations of flapping wings produce a camber effect that may modulate the production of lift and drag, and also alter aerodynamic control (Du and Sun, 2010; Walker et al., 2009; Zhao et al., 2009). Investigations of material properties and the behaviour of insect wings during cyclic loading are thus essential for our understanding of how elastic properties shape neuromuscular wing control, stroke kinematic patterns and aerodynamic force production to enhance aerodynamic efficiency in a flying insect. This study contributes to this research via an estimation of wing deformation power and demonstrates a source of mechanical energy loss due to spanwise wing bending in fly flight. The finding that the wing might recover only 80% bending energy during flapping motion, even in cases where aerodynamic power exceeds inertial, is not only relevant to understand the biology of flight but also for engineers who aim to improve aerodynamic performance and power expenditures of bio-inspired micro-air-vehicles.

APPENDIX

According to Walker (Walker, 2002) and correcting speed for wing translation and wing rotation, Magnus force (Fm) may be calculated by:
formula
(A1)
where ρ is the air density, c is the wing chord length, r is the distance from the wing hinge, is the angular velocity of wing translation, is the wing's morphological angle of attack with respect to the horizontal, and d is the differential. According to Ellington (Ellington, 1984b) and Lehmann and Pick (Lehmann and Pick, 2007), inertial forces due to translational acceleration of the wing in the horizontal () are given by:
formula
(A2)
where mw is the wing mass of each blade element and is the angular acceleration. Considering horizontal inertial forces generated by wing translation and wing rotation, we may modify the above equation to:
formula
(A3)
where l is the chordwise distance between the centre of wing mass of each blade and the rotational axis, and is the angular acceleration during wing rotation. Instantaneous inertial force in the vertical direction (Fi,v) depends on the wing's acceleration around the longitudinal axis and is thus:
formula
(A4)

Added mass force (Facc) is an acceleration reaction normal to the wing surface and may be expressed by the equation:

where is the blade velocity at 0.75c and normal to wing chord, and β is the added mass coefficient. Added mass reaction force of the fluid is already added to lift L and drag D estimates [eqns 6 and 7 in Walker (Walker, 2002); Fig. 5A,E]. Instantaneous aerodynamic force acting normal to the wing surface (Fa,n) equals the product between the vector sum of lift and drag and the wings' angular orientation (αg), given by:

formula
(A5)

Similar to the derivations above, we determined inertial forces acting normal to the wing surface (Fi,n) from the horizontal (Fi,h) and vertical (Fi,v) inertial components as:

formula
(A6)

The Magnus effect acts in the same direction of lift, thus the normal component of Magnus force (Fm,n) with respect to the wing is:

formula
(A7)

In Eqns A6, A7, A8, positive forces orientated normally to the wing surface pull on the dorsal wing surface and negative forces pull on the ventral surface.

formula
(A8)

Employing cantilever beam theory, we converted the measured spring constants of the fly wing into combined flexural stiffness using:

where k is the spring constant and EI is the mean flexural stiffness from the wing root to the measured point. Local flexural wing stiffness was derived by a method in which we considered the deflection of a blade at r as the difference between the local wing deflection and deflection caused by wing blades between r and the wing root. The deflection of a wing blade (δ; cantilever beam) when loaded at its end is given by:

formula
(A9)

where Fn is the force acting at the end of the wing blade. By contrast, the deflection of the blade end while loading a beam with length r>dr at r′ is given by the following expression:

formula
(A10)

where EI(r′) is the flexural stiffness at position rdr. By comparing the measured deflection at r using Eqn A10 with the deflection obtained from Eqn A11 at similar load, we obtain the blade deflection that is caused by the local flexural stiffness. Thus, we eventually derived the local EI from this deflection difference and using Eqn A10. We further used Eqn A11 to calculate wing tip deflection shown in Fig. 7. This measure is unequal to the sum of all local deflections because wing bending near the root causes larger wing tip deflection than bending near the wing tip because of the longer moment arm.

formula
(A11)

Local bending moment depends on the vector sum of aerodynamic force, added mass reaction force and inertial force normal to the wing surface multiplied by its moment arm. We ignored Magnus force because of previous suggestions (Walker, 2002). Thus local bending moment (M) is zero at the wing tip and maximum at the wing root, and equals the integral of all forces from distance r to the wing tip multiplied by their corresponding moment arms:

Total inertial power for wing motion (Pacc) was derived from the product between horizontal inertial force due to wing translation and rotation and the translational velocity of the blade, given by:

formula
(A12)

In contrast to Ellington's estimate of mean profile power requirements and instead of employing a drag coefficient based on Reynolds number, we estimated total aerodynamic costs (power, Paero) from translational, blade-wise wing velocity and our calculations of instantaneous horizontal drag (Eqn A3):

formula
(A13)
formula
(A14)

LIST OF SYMBOLS

     
  • c

    wing chord

  •  
  • mean wing chord

  •  
  • D

    drag parallel with local stream

  •  
  • E

    wing loading energy

  •  
  • EI

    flexural stiffness

  •  
  • EPE

    elastic potential energy

  •  
  • Fa

    total aerodynamic force

  •  
  • Facc

    added mass reaction force

  •  
  • Fa,h

    horizontal component of total aerodynamic force

  •  
  • Fa,n

    aerodynamic force normal to wing chord

  •  
  • Fa,v

    vertical component of total aerodynamic force

  •  
  • Fi

    total inertia

  •  
  • Fi,h

    inertia in the horizontal due to wing translation and rotation

  •  
  • inertia in the horizontal due to wing translation

  •  
  • inertia in the horizontal due to wing rotation

  •  
  • Fi,n

    total inertia normal to wing chord

  •  
  • Fi,v

    inertia in the vertical due to wing rotation

  •  
  • Fm

    Magnus force

  •  
  • Fm,n

    Magnus force normal to wing chord

  •  
  • Fn

    total force normal to wing chord (vector sum of

  •  
  • Fa,n,

    Fm,n and Fi,n)

  •  
  • h

    wing thickness

  •  
  • k

    combined wing spring constant from root to r

  •  
  • l

    distance of mw from rotational axis

  •  
  • L

    lift normal with local stream

  •  
  • mw

    wing mass of blade element

  •  
  • M

    local wing bending moment

  •  
  • n

    number of wing blade starting at wing root

  •  
  • Pacc

    total inertial power for wing flapping

  •  
  • Paero

    total aerodynamic power for wing flapping

  •  
  • r

    distance from wing root

  •  
  • r′

    blade position

  •  
  • R

    wing length

  •  
  • SED

    strain energy density

  •  
  • t

    time

  •  
  • v

    normal velocity of wing element computational

  •  
  • w

    width of wing blade

  •  
  • α

    geometrical angle of attack with respect to vertical

  •  
  • αg

    geometrical angle of attack with respect to horizontal

  •  
  • angular velocity of wing rotation

  •  
  • angular acceleration of wing rotation

  •  
  • β

    added mass coefficient

  •  
  • δ

    wing deflection

  •  
  • ρ

    air density

  •  
  • ε

    elasticity

  •  
  • ϕ

    wing stroke angle

  •  
  • angular velocity of wing translation

  •  
  • angular acceleration

FOOTNOTES

This work was funded by a scholarship of the Universiti Teknologi Malaysia to N.N. and grants LE905/9-3 and LE905/10-1 of the German Science Foundation to F.-O.L.

Acknowledgements

We would like to thank Ursula Seifert for editorial assistance.

REFERENCES

REFERENCES
Andersen
S. O.
,
Weis-Fogh
T.
(
1964
).
Resilin, a rubber-like protein in arthropod cuticle
.
Adv. Insect Physiol.
2
,
1
-
65
.
Balint
C. N.
,
Dickinson
M. H.
(
2004
).
Neuromuscular control of aerodynamic forces and moments in the blowfly, Calliphora vicina
.
J. Exp. Biol.
207
,
3813
-
3838
.
Bender
J. A.
,
Dickinson
M. H.
(
2006
).
A comparison of visual and haltere-mediated feedback in the control of body saccades in Drosophila melanogaster
.
J. Exp. Biol.
209
,
4597
-
4606
.
Biewener
A. A.
,
Roberts
T. J.
(
2000
).
Muscle and tendon contributions to force, work, and elastic energy savings: a comparative perspective
.
Exerc. Sports Sci. Rev.
28
,
99
-
107
.
Birch
J. M.
,
Dickinson
M. H.
(
2003
).
Spanwise flow and the attachment of the leading-edge vortex on insect wings
.
Nature
412
,
729
-
733
.
Bos
F. M.
,
Lentink
D.
,
Van Oudheusden
B. W.
,
Bijl
H.
(
2008
).
Influence of wing kinematics on aerodynamic performance in hovering insect flight
.
J. Fluid Mech.
5941
,
341
-
368
.
Casey
T. M.
,
Ellington
C. P.
(
1989
).
Energetics of insect flight
. In
Energy Transformations in Cells and Organisms
(ed.
Wieser
W.
,
Gnaiger
E.
), pp.
200
-
210
.
Stuttgart
:
Thieme
.
Combes
S. A.
,
Daniel
T. L.
(
2003a
).
Into thin air: contributions of aerodynamic and inertial-elastic forces to wing bending in the hawkmoth Manduca sexta
.
J. Exp. Biol.
206
,
2999
-
3006
.
Combes
S. A.
,
Daniel
T. L.
(
2003b
).
Flexural stiffness in insect wings I. Scaling and the influence of wing venation
.
J. Exp. Biol.
206
,
2979
-
2987
.
Combes
S. A.
,
Daniel
T. L.
(
2003c
).
Flexural stiffness in insect wings II. Spatial distribution and dynamic wing bending
.
J. Exp. Biol.
206
,
2989
-
2997
.
Dickinson
M. H.
,
Lighton
J. R. B.
(
1995
).
Muscle efficiency and elastic storage in the flight motor of Drosophila
.
Science
268
,
87
-
89
.
Dickinson
M. H.
,
Lehmann
F. O.
,
Sane
S. P.
(
1999
).
Wing rotation and the aerodynamic basis of insect flight
.
Science
284
,
1954
-
1960
.
Du
G.
,
Sun
M.
(
2010
).
Effects of wing deformation on aerodynamic forces in hovering hoverflies
.
J. Exp. Biol.
213
,
2273
-
2283
.
Ellington
C. P.
(
1984a
).
The aerodynamics of hovering insect flight. II. Morphological parameters
.
Philos. Trans. R. Soc. Lond. B
305
,
17
-
40
.
Ellington
C. P.
(
1984b
).
The aerodynamics of insect flight
.
III. Kinematics. Philos. Trans. R. Soc. Lond. B
305
,
41
-
78
.
Ellington
C. P.
(
1984c
).
The aerodynamics of insect flight. VI. Lift and power requirements
.
Philos. Trans. R. Soc. Lond. B
305
,
145
-
181
.
Ennos
A. R.
(
1988
).
The inertial cause of wing rotation in Diptera
.
J. Exp. Biol.
140
,
161
-
169
.
Ennos
A. R.
(
1989
).
Inertial and aerodynamic torques on the wings of Diptera in flight
.
J. Exp. Biol.
142
,
87
-
95
.
Fry
S. N.
,
Sayaman
R.
,
Dickinson
M. H.
(
2003
).
The aerodynamics of free-flight maneuvers in Drosophila
.
Science
300
,
495
-
498
.
Gorb
S. N.
(
1999
).
Serial elastic elements in the damselfly wing: mobile vein joints contain resilin
.
Naturwissenschaften
86
,
552
-
555
.
Haas
F.
,
Gorb
S. N.
,
Wootton
R. J.
(
2000a
).
Elastic joints in dermapteran hind wings: materials and wing folding
.
Arthropod Struct. Dev.
29
,
137
-
146
.
Haas
F.
,
Gorb
S. N.
,
Blickhan
R.
(
2000b
).
The function of resilin in beetle wings
.
Proc. R. Soc. Lond. B
267
,
1375
-
1381
.
Hedenström
A.
,
Ellington
C. P.
,
Wolf
T. J.
(
2001
).
Wing wear, aerodynamics and flight energetics in bumblebees (Bombus terrestris): an experimental study
.
Funct. Ecol.
15
,
417
-
422
.
Hepburn
H. R.
,
Chander
H. D.
(
1976
).
Material properties of arthropod cuticles: the arthrodial membranes
.
J. Comp. Physiol.
109
,
177
-
198
.
Ho
S.
,
Nassef
H.
,
Pornsinsirirak
N.
,
Tai
Y. C.
,
Ho
C. M.
(
2003
).
Unsteady aerodynamics and flow control for flapping wing flyers
.
Prog. Aerospace Sci.
39
,
635
-
681
.
Jensen
M.
,
Weis-Fogh
T.
(
1962
).
A rubber-like protein in insect cuticle
.
J. Exp. Biol.
37
,
889
-
907
.
Kesel
A. B.
,
Philippi
U.
,
Nachtigall
W.
(
1998
).
Biomechanical aspects of the insect wing: an analysis using the finite element method
.
Comput. Biol. Med.
281
,
423
-
437
.
Lehmann
F.-O.
(
2001
).
The efficiency of aerodynamic force production in Drosophila
.
J. Comp. Biochem. Physiol. A
131
,
77
-
88
.
Lehmann
F.-O.
,
Dickinson
M. H.
(
1997
).
The changes in power requirements and muscle efficiency during elevated force production in the fruit fly, Drosophila melanogaster
.
J. Exp. Biol.
200
,
1133
-
1143
.
Lehmann
F.-O.
,
Dickinson
M. H.
(
1998
).
The control of wing kinematics and flight forces in fruit flies (Drosophila spp.)
.
J. Exp. Biol.
201
,
385
-
401
.
Lehmann
F.-O.
,
Pick
S.
(
2007
).
The aerodynamic benefit of wing-wing interaction depends on stroke trajectory in flapping insect wings
.
J. Exp. Biol.
210
,
1362
-
1377
.
Lichtwark
G. A.
,
Barclay
C. J.
(
2010
).
The influence of tendon compliance on muscle power output and efficiency during cyclic contractions
.
J. Exp. Biol.
213
,
707
-
714
.
Liu
H.
,
Kawachi
K.
(
1998
).
A numerical study of insect flight
.
J. Comput. Phys.
146
,
124
-
156
.
Lv
S.
,
Dudek
D. M.
,
Cao
Y.
,
Balamurali
M. M.
,
Gosline
J.
,
Li
H.
(
2010
).
Designed biomaterials to mimic the mechanical properties of muscles
.
Nature
465
,
69
-
73
.
Maybury
W. J.
,
Lehmann
F.-O.
(
2004
).
The fluid dynamics of flight control by kinematic phase lag variation between two robotic insect wings
.
J. Exp. Biol.
207
,
4707
-
4726
.
Miller
L. A.
,
Peskin
C. S.
(
2005
).
A computational fluid dynamics of ‘clap and fling’ in the smallest insects
.
J. Exp. Biol.
208
,
195
-
212
.
Nachtigall
W.
(
1966
).
Die Kinematik der Schlagflügelbewegungen von Dipteren. Methodische und analytische Grundlagen zur Biophysik des Insektenflugs
.
Z. Vgl. Physiol.
52
,
155
-
211
.
Nachtigall
W.
,
Roth
W.
(
1983
).
Correlations between stationary measurable parameters of wing movement and aerodynamic force production the blowfly (Calliphora vicina R.-D.)
.
J. Comp. Physiol.
150
,
251
-
260
.
Neville
A. C.
(
1960
).
Aspect of the flight mechanics on anisopterous dragonflies
.
J. Exp. Biol.
37
,
631
-
656
.
Pfau
H. K.
(
1986
).
Untersuchungen zur Konstruktion, Funktion und Evolution des Flugapparates der Libellen (Insecta, Odonata)
.
Tijdschr. Entomol.
1293
,
35
-
123
.
Ramamurti
R.
,
Sandberg
W. C.
(
2002
).
A three-dimensional computational study of the aerodynamic mechanisms of insect flight
.
J. Exp. Biol.
205
,
1507
-
1518
.
Ramamurti
R.
,
Sandberg
W. C.
(
2007
).
A computational investigation of the three-dimensional unsteady aerodynamics of Drosophila hovering and manoeuvring
,
J. Exp. Biol.
210
,
881
-
896
.
Roberts
T. J.
,
Marsh
R. L.
,
Weyand
P. G.
,
Taylor
C. R.
(
1997
).
Muscular force in running turkeys: the economy of minimizing work
.
Science
275
,
1113
-
1115
.
Rüppell
G.
(
1989
).
Kinematic analysis of symmetrical flight manoeuvres of Odonata
.
J. Exp. Biol.
144
,
13
-
42
.
Sane
S. P.
,
Dickinson
M. H.
(
2001
).
The control of flight force by a flapping wing: lift and drag production
.
J. Exp. Biol.
204
,
2607
-
2626
.
Sane
S.
,
Dickinson
M. H.
(
2002
).
The aerodynamic effects of wing rotation and a revised quasi-steady model of flapping flight
.
J. Exp. Biol.
205
,
1087
-
1096
.
Shyy
W.
,
Liang
Y.
,
Tang
J.
,
Liu
H.
,
Trizila
O.
,
Stanford
B.
,
Bernal
L.
,
Cesnik
C.
,
Friedmann
P.
,
Ifju
P.
(
2008
).
Computational aerodynamics of low Reynolds number plunging, pitching and flexible wings
.
Acta Mech. Sinica
24
,
351
-
374
.
Smith
S. M.
(
1996
).
Biological control with trichogramma: advantage, successes and potential of their use
.
Annu. Rev. Entomol.
41
,
375
-
406
.
Steppan
S. J.
(
2000
).
Flexural stiffness patterns of butterfly wings (Papilionoidea)
.
J. Res. Lepidoptera
35
,
61
-
77
.
Sun
M.
,
Lan
S. L.
(
2004
).
A computational study of the aerodynamic forces and power requirements of dragonfly (Aeshna juncea) hovering
.
J. Exp. Biol.
207
,
1887
-
1901
.
Sun
M.
,
Tang
J.
(
2002
).
Unsteady aerodynamic force generation by a model fruit fly wing flapping motion
.
J. Exp. Biol.
205
,
55
-
70
.
Tu
M. S.
,
Dickinson
M. H.
(
1996
).
The control of wing kinematics by two steering muscles of the blowfly, Calliphora vicina
.
J. Comp. Physiol. A
178
,
813
-
830
.
Vanella
M.
,
Fitzgerald
T.
,
Preidikman
S.
,
Balaras
E.
,
Balachandran
B.
(
2009
).
Influence of flexibility on the aerodynamic performance of a hovering wing
.
J. Exp. Biol.
212
,
95
-
105
.
Walker
J. A.
(
2002
).
Rotational lift: something different or more of the same
.
J. Exp. Biol.
205
,
3783
-
3792
.
Walker
S. M.
,
Taylor
G. K.
,
Thomas
A. L. R.
(
2009
).
Deformable wing kinematics in the desert locust: how and why do camber, twist and topography vary through the stroke
.
J. R. Soc. Interface
6
,
735
-
747
.
Wang
H.
,
Zeng
L.
,
Liu
H.
,
Chunyong
Y.
(
2003
).
Measuring wing kinematics, flight trajectory and body attitude during forward flight and turning manoeuvres in dragonflies
.
J. Exp. Biol.
206
,
745
-
757
.
Wootton
R. J.
(
1991
).
The function morphology of the wing of Odonata
.
Adv. Odonatol
5
,
153
-
169
.
Wootton
R. J.
(
1992
).
Functional morphology of insect wings
.
Annu. Rev. Entomol.
37
,
113
-
140
.
Wootton
R. J.
(
1993
).
Leading edge section and asymmetric twisting in the wings of flying butterflies (Insecta, Papilionoidea)
.
J. Exp. Biol.
180
,
105
-
117
.
Wootton
R. J.
,
Evans
K. E.
,
Herbert
R.
,
Smith
C. W.
(
2000
).
The hind wing of the desert locust (Schistocerca Forskal). I. Functional morphology and mode of operation
.
J. Exp. Biol.
203
,
2921
-
2931
.
Young
J.
,
Walker
S. M.
,
Bomprey
R. J.
,
Taylor
G. K.
,
Thomas
A. L. R.
(
2009
).
Details of insect wing design and deformation enhance aerodynamic function and flight efficiency
.
Science
325
,
1549
-
1552
.
Zanker
J. M.
(
1990
).
The wing beat of Drosophila melanogaster
.
I. Kinematics. Philos. Trans. R. Soc. Lond. B
327
,
1
-
18
.
Zhao
L.
,
Huang
Q.
,
Deng
X.
,
Sane
P. S.
(
2009
).
Aerodynamic effects of flexibility in flapping wings
.
J. R. Soc.
7
,
485
-
497
.