Elastic deformation and energy loss of flapping fly wings

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 m² in lacewings, Hemerobius, to 10^–3 N m² 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 m², 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 m² whereas chordwise flexural stiffness ranged between 10^–5 and 10^–7 N m². 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.

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. x_r 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 mm², 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 mm³ (0.067), respectively.

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.

Wing elasticity (elastic efficiency, ε) is the ratio between unloading (E_U) and loading (E_L) 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).

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.

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.

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

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:

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.

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

The elastic force of blowfly wings linearly increased with increasing load applied chordwise until the displacement was constant (t₀, 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, R²=0.01, P=0.51, N=39; Fig. 4F).

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 m² 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 m², 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 m² (Fig. 4I).

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; F_i,h*=1.92 mN, F_i,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 (F_a,n) and inertia (F_i,n) normal to the wing surface amount to 474 and 235 μN, respectively, whereas added mass reaction force (166 μN) as part of F_a,n amounts to 35.0 and 70.6% of F_a,n and F_i,n, respectively. Similar values were obtained for the component of mean Magnus force acting normal to the wing surface (F_m,n=56.8 μN; 12.0% of mean F_a,n; 24.2% of mean F_i,n).

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 F_i,n=1.92 mN; maximum F_a,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 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, F_i,n=87.9±16.2 μN; time bin 0.4–0.5, F_i,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, F_a,n= 49.4 ±27.3 μN; time bin 0.4–0.5, F_a,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, F_m,n=–8.06±4.45 μN; time bin 0.4–0.5, F_m,n=8.65±4.78 μN; means ± s.d., N=9 wing blades, downstroke; Fig. 6B,F).

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

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

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×10³ J m^–3; 0.4–0.5 stroke cycle, 2.81±0.55×10³ 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).

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

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 m²) than our measurements (1.2×10^–7 N m²; 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.

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 A13–A14 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).

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.

Added mass force (F_acc) 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 (F_a,n) equals the product between the vector sum of lift and drag and the wings' angular orientation (α_g), given by:

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

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

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.

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:

where F_n 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:

where EI(r′) is the flexural stiffness at position r–dr. 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.

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 (P_acc) was derived from the product between horizontal inertial force due to wing translation and rotation and the translational velocity of the blade, given by:

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, P_aero) from translational, blade-wise wing velocity and our calculations of instantaneous horizontal drag (Eqn A3):

 
• c

  wing chord

 
• 

  mean wing chord

 
• D

  drag parallel with local stream

 
• E

  wing loading energy

 
• EI

  flexural stiffness

 
• EPE

  elastic potential energy

 
• F_a

  total aerodynamic force

 
• F_acc

  added mass reaction force

 
• F_a,h

  horizontal component of total aerodynamic force

 
• F_a,n

  aerodynamic force normal to wing chord

 
• F_a,v

  vertical component of total aerodynamic force

 
• F_i

  total inertia

 
• F_i,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

 
• F_i,n

  total inertia normal to wing chord

 
• F_i,v

  inertia in the vertical due to wing rotation

 
• F_m

  Magnus force

 
• F_m,n

  Magnus force normal to wing chord

 
• F_n

  total force normal to wing chord (vector sum of

 
• F_a,n,

  F_m,n and F_i,n)

 
• h

  wing thickness

 
• k

  combined wing spring constant from root to r

 
• l

  distance of m_w from rotational axis

 
• L

  lift normal with local stream

 
• m_w

  wing mass of blade element

 
• M

  local wing bending moment

 
• n

  number of wing blade starting at wing root

 
• P_acc

  total inertial power for wing flapping

 
• P_aero

  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

We would like to thank Ursula Seifert for editorial assistance.