A semi-empirical model of the aerodynamics of manoeuvring insect flight.

Blade element modelling provides a quick analytical method for estimating the aerodynamic forces produced during insect flight, but such models have yet to be tested rigorously using kinematic data recorded from free-flying insects. This is largely because of the paucity of detailed free-flight kinematic data, but also because analytical limitations in existing blade element models mean that they cannot incorporate the complex three-dimensional movements of the wings and body that occur during insect flight. Here, we present a blade element model with empirically fitted aerodynamic force coefficients that incorporates the full three-dimensional wing kinematics of manoeuvring Eristalis hoverflies, including torsional deformation of their wings. The two free parameters were fitted to a large free-flight dataset comprising N = 26 541 wingbeats, and the fitted model captured approximately 80% of the variation in the stroke-averaged forces in the sagittal plane. We tested the robustness of the model by subsampling the data, and found little variation in the parameter estimates across subsamples comprising 10% of the flight sequences. The simplicity and generality of the model that we present is such that it can be readily applied to kinematic datasets from other insects, and also used for the study of insect flight dynamics.

Introduction

The unsteady aerodynamics of insect flight have been the focus of considerable research, with new aerodynamic mechanisms still being discovered [1]. Dipteran flies have received particular attention, because their possession of only two functional wings reduces their kinematic complexity relative to four-winged insects, and avoids the aerodynamic complexity of tandem wing-wing interactions. Nevertheless, like most other insects, flies use high angles of attack and rapid wing rotation at stroke reversal, posing substantial challenges for aerodynamic modelling. Various modelling approaches have been used, each with their own advantages and disadvantages. More sophisticated techniques include the use of mechanical flappers [2-4] or computational fluid dynamics (CFD) (e.g. [5-8]). Both approaches determine the aerodynamic forces from a predefined set of wing kinematics, allowing the effect of specific kinematic parameters to be investigated experimentally [4,8,9], but recreating an insect’s wing kinematics in a mechanical or computational model is not straightforward, because the wings follow complex three-dimensional paths, and undergo substantial deformation through the wingbeat [10,11]. Wing deformation is difficult to replicate accurately in a mechanical flapper, and substantial user effort is required to generate a mesh capable of accommodating wing deformation in a computational model. In addition, the aerodynamics of the wings are affected by the body’s motions during flight manoeuvres, which are difficult to reproduce in a mechanical flapper and demanding to model computationally. Although the development of efficient algorithms and the increasing availability of low-cost clusters is making computational approaches ever-more practical, their application is still limited to quite small datasets.

A simpler approach is to use an analytical blade element model to estimate the aerodynamic forces, by splitting the wing into a series of narrow chordwise elements, each of which is modelled independently [12]. To the extent that the flow around a real wing is inherently three-dimensional and coupled to the wake [13], a blade element model cannot capture all of the details of the unsteady flow in the way that a computational model can. Nevertheless, as we demonstrate here, it is still possible to use this approach to make practically useful predictions of the forces by using empirically fitted force coefficients to summarize the complexities of the aerodynamics. Current blade element models comprise a quasi-steady component capturing how the pressure forces depend on the instantaneous velocity and angular velocity of the wing, and an unsteady component capturing how the pressure forces depend on the wing’s instantaneous acceleration [7,14,15], through the phenomenon of added mass. Other unsteady effects relating to the development of the flow are not captured by these models, because their coefficients are time-invariant and neglect the effects of wing–wake interactions from one half-stroke to the next [7]. Analytical blade element models therefore simplify the unsteady three-dimensional aerodynamics of flapping flight substantially, but can still do a surprisingly good job of approximating the forces produced by a flapping wing [7,14,15]. The key is to identify an appropriate analytical formulation describing how the aerodynamic forces on each blade element vary with the wing kinematics, and to model the wing kinematics with a sufficiently high degree of fidelity, which is the aim of this paper.

Previous studies have adopted several ad hoc approaches to modelling different aspects of the aerodynamics—especially the effects of wing rotation, which have been treated separately from the effects of wing translation in almost all previous models ([6,7,14,15]; but see [16]). In each case, these models are parametrized by a set of empirical force coefficients fitted to measurements made under different kinematic conditions using either a mechanical flapper [6,7,15] or a numerical model [14]. The resulting blade element models involve either one [6,7,15] or two [14] fitted parameters to describe the rotational lift and drag, together with separate expressions for the translational lift and drag coefficients as functions of the angle of attack, with two [14] to four [7,15] or five [6] fitted parameters. Consequently, the lift and drag are predicted from analytical expressions that are sometimes quite far removed from the underlying physics, and which involve from six [14] to nine [7,15] or even 11 [6] empirically fitted parameters. This brings an attendant risk of over-fitting, and as the identification and verification of these models has only been done using flat, rigid wings and simplified kinematics for the simplest case of equilibrium hovering flight [6,7,14,15], it is unknown how well they predict the aerodynamic forces and moments on real insects undergoing free-flight manoeuvres involving complex wing deformations. A particular concern with using such multi-parameter models of ad hoc form is that these may not generalize well to other flight morphologies, wing kinematics or flow conditions beyond those under which the data were collected. What is needed instead is a standard form developed from first principles—and therefore expected to generalize—that is capable of capturing the full complexity of the deforming wing kinematics.

In fact, as we show here, it is possible to predict approximately 80% of the variation in the stroke-averaged forces of free-flying hoverflies by using a physics-based model with just two free numerical parameters that are fitted empirically to the data. We achieve this by using linear least-squares modelling to fit the numerical coefficients of an unsteady blade element model, developed from first principles, to a free-flight dataset recording the deforming wing kinematics and stroke-averaged body accelerations for N = 26 541 wingbeats. Fitting this simple physics-based model to our free-flight dataset also allows us to interpret some of the more complicated empirical functions that have been fitted to model the aerodynamic force coefficients previously [3,6] and that have been co-opted into subsequent models [7,15]. Our new analytical blade element model takes full account of the three-dimensional motion of the wings and body, incorporating wing deformation in the form of a linear time-varying wing twist distribution, which is sufficient to capture most of the deformation that is present in Eristalis [11]. While our empirical data from hoverflies do not allow us to verify how accurately our model predicts the time history of the aerodynamic forces within a single wingbeat, the stroke-averaged forces are modelled closely. Given that the body dynamics of Eristalis are too slow to depend closely on the periodic forcing experienced through the wingbeat [17], our model’s ability to fit the stroke-averaged forces closely makes it well suited to use in future analyses of flight dynamics and control in this species. The simplicity and generality of our model is such that it can also be applied to kinematic datasets from other insects. Moreover, for the largest insects, in which the body dynamics operate on a similar timescale to the wingbeat [17], it should be possible to fit the aerodynamic force coefficients directly to the time-varying rather than stroke-averaged forces using the same method.

Methods

The overall aim of this paper is to develop an aerodynamic model with empirically fitted coefficients that predicts the stroke-averaged aerodynamic forces as a function of the instantaneous kinematic state of the wing through the stroke. The empirical measurements that we analyse here are those previously described in [18], but we begin by providing a brief summary of the experimental methods and kinematic reconstruction technique for context, before providing a detailed description of the aerodynamic modelling that forms the primary contribution of this paper.

Experimental methods

Adult Eristalis tenax and E. pertinax (Diptera: Syrphidae) were caught in Oxford and released singly inside a 1m diameter opaque acrylic sphere. Four synchronized high-speed video cameras (SA3, Photron Ltd, Bucks, UK) with 180 mm macro lenses (Sigma Imaging Ltd, Welwyn Garden City, UK) were used to record 768 × 640 pixel images at 3800 Hz. Bright back-illumination was provided by two synchronized 200W infrared pulsed lasers (HSI-5000, Oxford Lasers Ltd, Oxford, UK), each of which was routed through a split liquid light guide before being collimated by one of four large Fresnel lenses. The 805 nm wavelength of the laser illumination was far beyond the 600 nm upper limit of the visible spectrum for Eristalis [19-21], and a 20 μs pulse duration was used to eliminate motion blur and to prevent overheating of the insect. Ambient lighting was provided by an overhead LED light source. The cameras were self-calibrated using custom-written software in Matlab (The Mathworks Inc., Natick, MA, USA), to identify jointly optimal estimates of the camera parameters and the spatial coordinates of points on a two-dimensional calibration grid held in a range of positions and orientations [22]. We captured between 10 and 50 separate recordings of each hoverfly, before anaesthetizing it with CO2 at the end of the experiment, and weighing it using a microbalance with 0.1 μg readability (UMX2, Mettler Toledo Ltd, Leicester, UK). We measured the temperature (27.5 ± 0.9°C) and relative humidity (33.4 ± 4.0%) of the air in each trial (mean ± s.d. for 36 insects) and used these to determine the air density (ρ). We analysed all of the flight sequences in which both wings were visible to all four cameras for ≥2 wingbeats, giving a total sample of N = 26 541 wingbeat pairs from 854 flight sequences representing 36 hoverflies. A fully automated shape-carving procedure was used to label voxels contained within the minimum convex hulls of the insect’s body and its two wing outlines, respectively [18].

Flight kinematics modelling

Throughout the paper, we use boldface symbols to represent vector quantities, and use t to represent continuous time. We defined the body kinematics using a right-handed, rotating, body-fixed axis system {x, y, z} with its origin at the centre of volume of the body voxels, its x-axis directed anteriorly along the major axis of the body voxels, and its y-axis pointing rightward parallel to the line connecting the wing bases (figure 1). We measured the position of the body axes, (t), in a non-rotating laboratory coordinate system {X, Y, Z}, in which the Z-axis was vertical, and smoothed our measurements using a quintic spline fitted in B-form in Matlab [23]. This method fits each element of the position vector as the smoothest piecewise polynomial function of time that falls within a given tolerance of the data, defined as the sum of the squared distance of the function from the data over all sample points. For transparency, we selected a tolerance that was equal to the sum of the squared variation that would have been removed by forward–backward filtering the data using a third-order Butterworth filter with a 100 Hz cut-off frequency (−3 dB) chosen to fall well below the insect’s 188 ± 14 Hz wingbeat frequency (mean ± s.d.). We then double-differentiated each quintic spline analytically in Matlab to estimate the insect’s instantaneous linear acceleration with respect to the laboratory coordinate system, , resolved in the laboratory axes.

Figure 1

Kinematic definitions. (a) The body axes {x, y, z} are defined as a right-handed axis system (grey arrows) with its origin at the centre of volume of the insect’s body (grey silhouette), its x-axis directed anteriorly along the major axis of the body voxels, and its y-axis pointing to the insect’s right, parallel to the line connecting the wing bases. The tip kinematics of the right wing are defined by its spherical coordinates in another right-handed axis system {x, y, z} (black arrows) fixed parallel to the body axes, but originating at the wing base. The tip kinematics of the left wing are defined by its spherical coordinates in an equivalent left-handed axis system {x, y, z} (not shown). In each case, the line from wing base to wing tip (red dot) defines the spanwise rotational axis of the wing (red arrow): the azimuth of this spanwise axis defines the stroke angle of the wing (ϕ), and its elevation defines the deviation angle (θ). The blue shaded area shows a single chordwise blade element with the position its three-quarter chord point marked by a green dot. (b) The pitch angle ω of a blade element is defined having first rotated the wing’s measured outline through its stroke angle −ϕ about the z-axis, then through its deviation angle −θ about the x-axis, so as to bring the line from wing base to wing tip into alignment with the y-axis. The pitch angle ω was then defined as the angle from the xy-plane to the anatomical ventral side of the rotated chord perpendicular to the y-axis. (c) The speed of a blade element (U) is measured at its three-quarter chord point (green dot), and is the hypotenuse of the components of the velocity vector directed parallel to the blade element chord () and normal to the blade element surface (U⊥). The aerodynamic angle of attack (α) is defined as shown by the arctangent of these two vector components.

We described the kinematics of the right wing in a right-handed axis system {x, y, z} parallel to the rotating body axes {x, y, z}, but with its origin at the wing base. We used the spherical coordinates of the wing tip to define the stroke angle (ϕ) and deviation angle (θ) of the wing (figure 1). We defined the local pitch angle ω(r) of the wing at radial coordinate r (figure 2) by rotating the wing’s measured outline through the angle −ϕ about z, then through the angle −θ about x, so as to bring its spanwise axis into alignment with y. The local pitch angle ω(r) of the wing was then defined as the angle from the xy-plane to the anatomical ventral side of the chord, perpendicular to the y-axis at radial coordinate r (figures 1b and 2). We summarized the spanwise twist as a linear function of distance along the wing by regressing ω(r) on r, modelling the local pitch angle as:where ω0 is the pitch angle offset, and ω is the linear twist gradient. The kinematics of the left wing were defined independently using an equivalent left-handed axis system {x, y, z}. Variable wing camber might also have been present [11], but measuring this requires data of higher order than can be obtained using a voxel carving method to identify the wing outlines. Were wing camber to be measured directly in a future study, it would in principle be possible to incorporate its effects in the blade element model below by treating camber as augmenting the aerodynamic angle of attack defined with respect to the angle of the chord line measured between the leading and trailing edge.

Figure 2

Morphological definitions and measurements. (a) The wing is split into 20 equally spaced chordwise blade elements, each of width Δr. The ith blade element is located at radial distance r(i) from the wing base, and has chord length c(i). The distance of the three-quarter chord point (green circle) from the spanwise rotation axis of the wing (red arrow) is denoted c3/4(i), and is used in calculating the quasi-steady lift and drag forces (see electronic supplementary material, Methods). Likewise, the distance of the half-chord point (black circle) from the spanwise rotation axis is denoted c1/2(i), and is used in calculating the added mass force (see electronic supplementary material, Methods). (b) The measured positions of the three-quarter chord points (green circles) and half-chord points (black circles) are shown for each of the 20 blade elements.

We smoothed our measurements of θ, ϕ, ω and ω0 for each wing using the same quintic spline method as we used to smooth the body kinematics, this time setting the tolerance to give a similar degree of smoothing to a digital Butterworth filter with a cutoff frequency of 500 Hz for the wing tip kinematics θ and ϕ and 800 Hz for the wing twist kinematics ω and ω0 (electronic supplementary material, figure S1). We then evaluated the first and second derivatives of these fitted spline functions analytically. The smoothed kinematic data were next split into discrete wingbeats by identifying the time at which the mean angular speed of the two wing tips reached a minimum at the end of each half stroke. Finally, for consistency of sampling between discrete wingbeats of variable period, we used cubic interpolation to resample the smoothed wing and body kinematic data and their time derivatives at 100 evenly spaced time steps through each wingbeat, beginning at the start of the downstroke. This means that the kinematic measurements were upsampled by approximately a factor of 5 prior to further analysis, which (i) ensures that each wingbeat starts and ends at exactly the same phase; (ii) allows the data to be stored in an efficient matrix form; and (iii) standardizes the basis on which the wingbeat averaged forces are estimated.

By the end of this process, each wingbeat is represented by 4200 datapoints, comprising the 6 degrees of freedom of the body and four primary kinematic variables of each wing, each sampled together with their first and second derivatives 100 times per wingbeat. In principle, these data are already in a form suitable for the subsequent blade element analysis. However, in the interests of compressing the data into a compact functional form suitable for sharing, and noting the different characteristic timescales on which the different kinematic components vary, we projected the data for each wingbeat into a set of harmonic basis functions [24] comprising a truncated Fourier series plus cubic polynomial in time from the start of the wingbeat, with harmonic content to first order for the body kinematics, fourth order for the wing tip kinematics and sixth order for the wing twist kinematics (see electronic supplementary material, Methods). This compression reduces the dimension of the data by almost a factor of 30, while preserving >99.99% of the measured variation in the pose of the insect, as characterized by the 6 degrees of freedom of its body and four primary kinematic variables of each of its wings. These harmonic representations of the data are shared as electronic supplementary material, Data S1, so to ensure the repeatability of our analysis and to enable its validation using numerical techniques possibly requiring finer time steps, we obtained the 4200 datapoints that we use in the blade element model for each wingbeat by evaluating these harmonic fits rather than the spline fits on which they are modelled. This step makes a negligible difference to the numerical values of the datapoints used in the analysis (electronic supplementary material, figure S1), and hence to the results of the analysis, but it aligns the present work more closely to the approaches that we have developed elsewhere for analysing the dominant kinematic couplings involved in insect flight control using harmonic functional principal components analysis [24].

Standard hovering wingbeat

We defined a standard wingbeat for use in model validation by taking the mean through time of a sample comprising the 1% of all N = 26 541 wingbeat pairs that most closely met the criteria for hovering flight. We selected these as the 265 wingbeats with the lowest flight speed from within the set of wingbeats representing near-equilibrium flight, which we defined as flight where the magnitude of the body’s acceleration was less than 0.5 m s−2 in both the vertical and the horizontal (i.e. such that the vertical aerodynamic force would have been within 5% of supporting body weight). The 265 wingbeats that we used in this averaging came from 51 different flight recordings and 19 different individuals, and the standard hovering wingbeat that they define should therefore be representative of equilibrium hovering flight (figure 5). The cubic polynomial terms of the harmonic representation of the wing kinematics (see above) were negligible for the standard hovering wingbeat, and we therefore set them to zero so as to make the standard hovering wingbeat kinematics strictly periodic. The body was assumed to be stationary, with a constant pitch angle of 43.6° between the long axis of the body and the horizontal, defined as the mean over the 265 wingbeats, and with the body roll angle set to zero as appropriate for symmetric hovering.

Flight dynamics modelling

The overall goal of this paper is to use aerodynamic modelling to relate free-flight measurements of body kinematics to the wing kinematics that produce them. Although we take full account of the body’s rotational and translational motion in defining the motion of the wings relative to the air, our modelling of the resulting aerodynamic forces only considers their effects on the translational motion of the centre of mass. That is to say, we do not attempt to model the rotational dynamics of the body, which is a more complex problem requiring knowledge of the insect’s inertia tensor and the chordwise position of the centre of pressure. Subject to making such further assumptions, the rotational dynamics can be modelled subsequently using the same blade element model if required.

The only ways in which a fluid can impart force to a solid surface are through pressure forces acting normal to the surface, and friction forces acting tangential to it. We may therefore use Newton’s Second Law to write the equations of translational motion for a free-flying insect aswhere m is the insect’s mass, is the acceleration of the insect’s centre of mass with respect to the laboratory coordinate system, is gravitational acceleration, is the total pressure force, and is the total friction force. Note that as equation (2.2) contains only the external forces acting at the insect’s centre of mass, it does not show the inertial forces that act in reaction to flapping at the wing hinge. These are internal forces that cannot produce any acceleration of the insect’s centre of mass, which is therefore a different situation to that which is encountered when measuring the internal forces at the wing hinge on a bench-mounted mechanical flapper. Moreover, although the wings’ motion causes the anatomical position of the centre of mass to fluctuate, which can cause the body to oscillate at wingbeat frequency in some large insects [25], the wings are several orders of magnitude lighter than the body in Eristalis. Their motion therefore has a negligible effect on the motion of the centre of volume of the body voxels, which we thereby equate with the motion of the insect’s centre of mass.

Since the wingbeat period of a hoverfly is much shorter than any characteristic timescale of its body’s dynamics [17], we may reasonably model its body dynamics using the stroke-averaged versions of these variables, which we denote using overbar notation with n as wingbeat number:This stroke-averaging is beneficial in removing the noise associated with estimating the double derivative of position from high-speed videogrammetric data. Nevertheless, in future work with larger insects whose body dynamics operate on a similar timescale to the wingbeat [17], it might be more appropriate to retain the original form of equation (2.2). The left-hand side of equation (2.3) is a direct empirical estimate of the stroke-averaged aerodynamic forces, which we refer to hereon as the ‘measured’ aerodynamic force, to distinguish this from the indirect estimates of the time-varying aerodynamic force that we make later using our blade element model. The stroke-averaged expression in equation (2.3) forms the basis of the empirical estimation of the aerodynamic forces in this paper.

Aerodynamic modelling

We begin this section by analysing the scaling of the aerodynamic forces, which we use to define the theoretical form of the kinematic predictor variables that we use to fit the aerodynamic coefficients of the quasi-steady blade element model. Although the key concepts are covered in aerodynamics texts such as Katz & Plotkin [26] and in the seminal work on insect flight by Ellington [12], they are usually only discussed in the context of detailed aerodynamic models that aim to determine analytically the same force coefficients as we aim to estimate empirically here, under restrictive assumptions that will not usually be satisfied in insect flight. Hence, following the approach advocated by [27], and noting that there is no analytical theory which covers the full range of flow conditions experienced in insect flight, we here develop the scaling of the aerodynamic forces from first principles, with the goal of making the assumptions of our blade element model fully transparent. We then introduce the theoretical form of the equations that we use to model the unsteady forces, before describing how we fit the empirical aerodynamic force coefficients to the measured flight data. Although some elements of our model are shared with previous work, this approach results in an aerodynamic model with a different mathematical form to the others that have been used previously to analyse insect flight [3,6,7,12,14,15].

Scaling of the quasi-steady pressure force

The pressure force in equation (2.2) will dominate the friction force at the Reynolds numbers of order 103 that characterize flight in Eristalis. It follows that viscous shear in the thin boundary layer forming the inner region of the flow cannot be directly responsible for setting the fluid into motion in the outer region of the flow. This outer flow must instead be driven by pressure gradients resulting from the wing’s motion and reflecting the constraint that fluid cannot penetrate its surface. This constraint implies that at every point on the wing’s surface, the surface-normal component of the local flow velocity (Q⊥) must be identical to the surface-normal component of the local surface velocity (V⊥), each measured with respect to an inertial frame of reference:where V is the speed of a given surface point, and where β is the angle of attack measured from the velocity of this point to the plane tangent to the surface there. This boundary condition is not sufficient to uniquely determine the tangential flow [26], but the pressure gradients that result from its imposition are nevertheless causative in generating the tangential flow. For example, the line integral of the tangential flow around the aerofoil is called the circulation (Γ), and if it is assumed that viscous effects cause the flow to depart smoothly from the sharp trailing edge and to have finite velocity there, which is called the Kutta condition, then the circulation that is needed to maintain this condition at steady state may be expected to scale with the normal flow that would be present at the trailing edge in the absence of any circulation. A similar conclusion holds in relation to a generalized form of the Kutta condition that can be applied at a sharp leading edge, wherein the normal velocity of the flow determines both the critical angle of attack at which the flow begins to separate, and the strength of the resulting leading-edge vortex [28]. To summarize, a moving wing perturbs the flow because of its impermeability to the fluid, and the form of the resulting kinematic boundary condition (equation 2.4) immediately suggests that the strength of this perturbation should scale with the speed of the surface and the sine of its angle of attack.

To specialize this still quite-general conclusion, we note that for a thin flat plate undergoing steady translational motion through still air, the local speed V and angle of attack β is the same everywhere. Given that the circulation is defined as the line integral of the tangential flow around the chord of the plate, we may therefore expect it to scale as where U and α are used to represent the overall speed and angle of attack of the aerofoil, and where c is its chord length. This conclusion accords with the results of classical thin aerofoil theory, which predicts that for a thin flat plate undergoing steady translational motion with the Kutta condition met at the trailing edge [26]. Of course, if the aerofoil is also rotating steadily, then the local angle of attack β will vary linearly along the chord, so it is reasonable to ask whether the quasi-steady effects of translation and rotation may in fact be captured by equating the overall angle of attack α to the local angle of attack β measured at some unique point on the chord. Indeed, in the limit of small amplitude harmonic oscillations, classical unsteady thin aerofoil theory predicts that the circulation is independent of the angular velocity of the aerofoil in the special case that the angle of attack is measured at the three-quarter chord point [26]; and is thereby proportional to the normal flow velocity there [28]. Here, we test empirically whether this conclusion also holds for our data.

The scaling of the circulatory pressure force is given by the Kutta–Joukowski theorem for two-dimensional inviscid flows as , where ρ is the fluid density [26]. The same theorem also predicts that this pressure force will act perpendicular to the flow—not perpendicular to the chord. This is unintuitive for a thin flat plate generating a pressure difference across its surfaces, but reflects the fact that if the flow remains attached then the acceleration of the fluid around the leading edge will be associated with a suction force parallel to the chord. In practice, this leading-edge suction is expected to be lost if the flow separates at the leading edge, as is typical of insect flight, and it is then more reasonable to assume that the pressure force will act normal to the chord. Furthermore, on a three-dimensional wing of low aspect ratio, the flow induced by the wake tilts the aerodynamic force vector back considerably, bringing its line of action more nearly perpendicular to the chord. Putting all of these considerations together, we therefore propose modelling the quasi-steady pressure force (Pqs) on a blade element of small width Δr using the scaling:where the unknown constant of proportionality and direction of action of the force remain to be determined empirically. Decomposing this quasi-steady pressure force into a lift component L = Pqscos α, and a drag component D = Pqssin α then yields the following scalings for the pressure lift and drag:andwhere the constant of proportionality is assumed to be identical in both expressions, and where by definition the drag acts in the direction of the flow and the lift acts perpendicular to the flow in the direction of increasing angle of attack. We emphasize that these equations are only valid where the angle of attack α is measured at a unique reference point on the chord, provisionally assumed to be located at three-quarters of the distance from the leading edge to the trailing edge, where α is defined as the angle measured between the chordline and the velocity of that point on the chord.

Quasi-steady force coefficients

By convention, the non-dimensional lift coefficient (C) and drag coefficient (C) of a blade element are defined asandCombining these identities with the scalings in equations (2.6) and (2.7), we expect that the parts of the lift and drag coefficients due to the quasi-steady pressure force should vary asandwhere is an unknown constant that remains to be determined. This parameter may be interpreted as the derivative of C with respect to the angle of attack α at α = 0, where C is the force coefficient associated with the total quasi-steady pressure force Pqs, non-dimensionalized similarly to C and C in equations (2.8) and (2.9).

Equation (2.11) implies that the drag force should vanish at α = 0, which reflects our neglect of the friction force to this point. However, whereas the pressure force is expected to dominate the friction force at most angles of attack, friction cannot be dismissed entirely at very low angles of attack, when the pressure force will also be small (equation (2.5)). Because friction acts tangential to the surface, it will in principle modify the total lift and drag coefficients asandwhere the notation C(α) reflects the fact that the dimensionless friction force coefficient C will in general depend on the angle of attack α. However, on the basis that the friction forces are expected to be negligible in comparison to the pressure forces except as α → 0, when C(α) sin α → 0 and C(α) cos α → C(α), it is reasonable to model the lift and drag coefficients empirically asandwhere the tilde notation indicates an empirical parameter estimate, and where the drag coefficient offset is equal to the drag coefficient at α = 0 such that . We note in passing that this formulation is mathematically equivalent to assuming that the friction coefficient varies as , which appropriately predicts that the net tangential force on an infinitely thin flat plate will vanish at α = 90°. With this assumption, it follows that , and it can be seen by inspection of equations (2.12) and (2.14) that , where the Pythagorean identity sin2 α + cos2 α = 1 is used to prove the same in equations (2.13) and (2.15). The estimated lift and drag coefficients are the same on either basis, and the only distinction to be drawn here is in whether is interpreted as an empirical estimate of or of .

This approach neglects possible variation of the friction drag coefficient with Reynolds number through the stroke, but this is unlikely to represent much of a limitation in practice, given that the drag coefficient offset is determined empirically, and lumps the friction drag together with any pressure drag that may happen to be present at α = 0 owing to the effects of wing corrugations, leading-edge thickness, etc. Such effects are difficult to model from first principles, being highly dependent on the detailed wing structure.

Form of the unsteady forces

Differentiating the boundary condition in equation (2.4) with respect to time shows that the acceleration of any point on the wing normal to its surface must be accompanied by an equal acceleration of the fluid at that point. Pressure will communicate this motion to the rest of the fluid at the speed of sound, decaying quickly away from the wing, and effectively superposing instantaneously with the existing flow. The reaction to this acceleration of the fluid is called the added mass force, because it has the effect of increasing the inertia experienced by a wing accelerating through a fluid. It is clear from the form of the boundary condition, however, that in contrast to the force that is required to accelerate the mass of a solid body, the added mass force depends on how the body’s acceleration is directed in relation to its geometry. Nevertheless, by resolving the body’s acceleration relative to the inertial frame in body axes [28], it becomes possible to define an added mass tensor that depends only on the shape and size of the body. The added mass forces involved in insect flight are especially significant at stroke reversal [29], when the magnitude of the wing’s acceleration is greatest, and when the wing’s orientation is such that a significant part of its acceleration is directed normal to its chord. It is important to emphasize that the added mass force is simply an unsteady pressure force arising from acceleration of the fluid, and is not connected with any identifiable mass of fluid actually being carried along with the body.

Because the instantaneous effect of the wing’s acceleration superposes with the existing flow, inviscid flow theory can be used to model the added mass forces acting in a real fluid. Indeed, it has been established theoretically [30], validated numerically [29,30], and demonstrated experimentally [31,32], that the same added mass force acts in both viscous and inviscid flow. Classical inviscid flow theory [33,34] predicts the surface-normal component of the added mass force on a thin flat plate of chord c and width Δr aswhere denotes the scalar projection of the acceleration of the plate’s half-chord point onto the unit vector normal to the plate’s surface. We emphasize that the validity of this equation rests on the assumption that the surface-normal acceleration is measured relative to an inertial frame of reference at the half-chord point. Equation (2.16) has been widely used to model the added mass forces on insect wings [7,12,14,15], but it is perhaps less well known that the same inviscid flow theory [33,34] also predicts a surface-parallel component of the added mass force, acting along the chord aswhere U′⊥ denotes the surface-normal velocity of the half-chord point, and Ω denotes the angular velocity of the plate along its spanwise axis [29,30]. It can be seen from the form of equation (2.17) that this surface-parallel component of the added mass force is solely a reaction to the centripetal part of the acceleration of the fluid parallel to the wing’s chord (i.e. there is no added mass force associated with tangential acceleration of a thin flat plate parallel to its chord). It is therefore directly attributable to the surface-normal motion of the plate, as expected from the boundary condition in equation (2.4). For a flat plate of finite width Δr, the added mass force will also include a surface-parallel component acting along the span, with the same form as equation (2.17) after switching the spanwise and chordwise axes (e.g. [33]).

For the general planar motion of a rigid aerofoil described by equations (2.17) and (2.16), the combined contributions of A⊥ and will integrate to zero over a periodic stroke cycle when resolved in inertial axes [34]. It is important to note, however, that the surface-normal component alone need not integrate to zero when resolved in inertial axes. This is the situation that arises if the wing rotates and the surface-parallel components of the added mass force are neglected—as they are in most blade element models [7,14,15]. This may not be unreasonable, because the flow evolves such that the circulatory force cancels the surface-parallel component of the added mass force when the flow separates at the leading edge, thereby causing the total pressure force on a rotating and translating flat plate to act approximately normal to its surface [30]. Flow separation is a viscous phenomenon, so it seems plausible that this cancellation of the surface-parallel component of the inviscid added mass force might be the origin of the viscous ‘centripetal acceleration reaction’ force reported by Zhang et al. [35], given that the contribution of the remaining surface-normal component need not integrate to zero. We therefore follow the usual convention of modelling only the surface-normal component of the added mass force here. This means that the added mass forces that we calculate do not integrate to zero over a single wingbeat period—though there would be no reason to expect them to do so anyway during stroke cycles that are not strictly periodic, such as those used during unsteady manoeuvres. Finally, we note that the classical inviscid flow theory cited above refers to single rigid bodies [33,34], and it is well established mathematically that deformable bodies [36] and articulated bodies [37] are capable of generating periodic locomotion through perfect fluids, solely by means of inviscid added mass effects [38].

Blade element modelling

The final step in the aerodynamic modelling is to use the equations in §2.5.2 to assemble a set of kinematic predictors for the stroke-averaged quasi-steady aerodynamic forces. Given our measurements of the wing and body kinematics, these results can then be used together with equations (2.2) and (2.16) to formulate a set of linear equations that can be solved in a least squares sense for the unknown aerodynamic force coefficients and , which are assumed to be the same for all blade elements. Practically speaking, we split each wing into 20 evenly spaced blade elements, of width Δr and chord length c(i), where i ∈ 1 … 20 denotes the blade element number (figure 2), the aerodynamic contributions of which we then summed over all 20 blade elements and all 100 sample points for both wings.

Making use of the equations in §2.5.2, we modelled the instantaneous lift (), drag () and added mass () forces acting on the ith blade element at time t as

andin which , 1⊥, 1⊥ are unit vectors defined as follows. The unit vector 1⊥ is directed mutually perpendicular to the velocity of the blade element and its span, and points in the direction of increasing angle of attack, thereby defining the direction in which lift acts. The unit vector is directed parallel to the velocity of the blade element, and points in the direction of the relative flow, thereby defining the direction in which drag acts. In each case, the blade-element velocity is defined at the three-quarter chord point. The unit vector 1⊥ is the blade-element surface normal, signed positive in the direction of increasing angle of attack, and is related to the other two unit vectors by the identity . The aerodynamic speed, U, and angle of attack, α, of each blade element are each measured at the three-quarter chord point working back from the leading edge, as detailed in the electronic supplementary material, Methods, such that equations (2.18) and (2.19) account implicitly for the effect of wing rotation on the quasi-steady pressure force (see also §2.5.1). The surface-normal acceleration is defined at the half-chord point for the purposes of determining the added mass force, and is signed positive in the direction of 1⊥ (see electronic supplementary material, Methods). All of these quantities measure the blade element kinematics with respect to the laboratory coordinate system, and therefore incorporate all of the effects of the body’s translational and rotational motion during forward flight and manoeuvring.

Our kinematic measurements record the motion of each wing at τ ∈ 1 … 100 discrete time points through each wingbeat. The contribution of each wing to the total stroke-averaged aerodynamic force may therefore be written as where we have made use of the identity to eliminate two of the trigonometric terms when combining equations (2.18) and (2.19). By an obvious use of notation for the summations, we will abbreviate the right-hand side of equation (2.21) as . The total stroke-averaged aerodynamic force on a given wingbeat may therefore be modelled aswhere the L and R superscripts denote summations for the left and right wings, respectively.

Model fitting

Combining equation (2.22) with equation (2.3) for the measured stroke-averaged aerodynamic force, we may writewhere and so on, and where is a residual error term accounting for the difference between the measurements and the model in the lab axes in which the forces are resolved. Equation (2.23) is linear in the unknown force coefficients and , so can be solved using linear regression to provide parameter estimates and that minimize the error sum of squares over all N = 26 541 wingbeats. We initially solved equation (2.23) with the aerodynamic speed U and angle of attack α defined at the three-quarter chord point as explained above, which yields a unique solution for and . As a direct check on the effect of this assumption of the model, we then tried varying the chordwise position at which U and α were defined, which yields a family of solutions for and . Finally, we verified the importance of including body motion and wing twist in the model by re-computing the kinematics without accounting for body motion, and with the wing pitch angle set uniformly at its value mid-span, before solving again for the unknown force coefficients and (table 1).

Table 1

Comparison of the full aerodynamic model with several alternative models, showing the effect of various simplifications on the estimated force coefficients and , and on the mean squared error (MSE) in the predicted forces. The penultimate column reports the chordwise reference point at which the kinematics are defined for the purposes of calculating the quasi-steady lift and drag on each blade element; the added mass forces are always calculated using the acceleration measured at 50% of chord. This chordwise reference point is optimized in the version of the full model shown in the last row of the table, so as to minimize the MSE subject to the necessary physical constraint that . The last column gives the total number of free parameters estimated by the optimization procedure. The MSE is averaged over all N wingbeats and all three axes such that , where contains the residual error in each axis for the nth wingbeat. The estimated force coefficients are reported ±1 s.d., where s.d. is the standard deviation of their values computed over 105 random subsamples of the data, each comprising 10% of the recorded flight sequences.

model	C^_P_α	C^_D_0	MSE (mN)2	% chord	parameters
full model	2.70 ±0.04	0.15 ± 0.07	0.0265	75	2
no drag offset	2.73 ± 0.04	—	0.0267	75	1
no drag offset or added mass	3.12 ± 0.03	—	0.0280	75	1
flat plate wing	1.69 ± 0.04	0.48 ± 0.06	0.0344	75	2
no body motion	2.65 ± 0.04	0.86 ± 0.19	0.0356	75	2
full model, optimized % chord	2.58 ± 0.04	0	0.0262	87	3

Model validation

In the light of the very large sample, and because the regression model does not take account of autocorrelation in the stroke-averaged forces from one wingbeat to the next, we do not report 95% confidence intervals for our parameter estimates. Instead, we tested the robustness of the analysis using subsampling, repeating the regression modelling 105 times on random subsamples of the data each containing only 10% of the flight sequences (i.e. 85 out of the 854 flight recorded sequences). This subsampling analysis allows us to assess the variance in our parameter estimates resulting from variation between individuals and flight sequences, and does so at a sample size that is more realistic for future studies than the very large sample used here (i.e. order 102 rather than order 103 flight sequences). Finally, using our estimates of and with the kinematics defined at the three-quarter chord point, we tested how the number of blade elements and number of time steps affected the aerodynamic forces predicted for the standard hovering wingbeat. The predicted mean absolute aerodynamic force changed by only 0.30% when increasing the number of blade elements to 1000 (electronic supplementary material, figure S2a,b), and by only 0.14% when increasing the number of time steps to 10 000 (electronic supplementary material, figure S2c,d), so we conclude that the default discretization using 20 blade elements and 100 time steps is more than sufficient.

Results

Body dynamics

Our free-flight dataset captures a wide variety of behaviours, including forward flight, hovering, ascent, descent and saccadic manoeuvres. A typical flight recording (electronic supplementary material, Video S1) includes brief periods of slow forward flight, punctuated by fast body saccades. Although we cannot claim to have captured the entire flight envelope, these data therefore cover a large part of the behavioural repertoire of Eristalis, including many manoeuvres typical of free-flight [39-42]. This range of behaviour is reflected in the variability of the measured aerodynamic forces acting along the x-axis (0.97 ± 0.43 mN) and z-axis (−1.11 ± 0.4 mN) of the body (mean ± s.d.; figure 3a,c). The forces measured along the y-axis were much less variable (0.00 ± 0.11 mN; figure 3b), consistent with the orthodoxy that comparatively little lateral aerodynamic force is produced during manoeuvres [43]. Nevertheless, the double derivatives of body position vary to a similar extent in all three body axes (figure 3d–f), indicating that the asymmetry of force production in the x- and z-axes is explained by the need to overcome the body’s acceleration due to gravity, which is small in the y-axis except during highly banked turns. It follows that after addressing the requirement for weight support, the hoverflies were actually comparably manoeuvrable in all three body axes.

Figure 3

Histograms of the measured stroke-averaged aerodynamic force and measured acceleration of the insect’s body in an inertial frame of reference. (a–c) Measured stroke-averaged aerodynamic force, resolved in the insect’s body axes. (d–f) Measured stroke-averaged acceleration, , resolved in the insect’s body axes. Note that the distribution of the acceleration is similar across all three body axes, but that the need to provide weight support in addition to manoeuvring force means that the resultant forces are principally distributed in the x- and z-axes of the body.

Wing kinematics

Wingbeat frequency and stroke amplitude vary greatly over the dataset, but their variability owes more to variation between individuals than within (figure 4), and the time-history of the wing kinematics is actually quite stereotyped over the whole dataset (figure 5b). Because the wing tip trajectory is inclined at approximately 45° to the body, the stroke angle ϕ and deviation angle θ always have similar oscillation amplitudes, and both vary approximately sinusoidally through the wingbeat (figure 5a). The wing pitch angle ω at mid-span varies symmetrically on the upstroke and downstroke, showing a slight recoil at the start of each half-stroke. The aerodynamic angle of attack α has a similar time history on both the upstroke and the downstroke, changing rapidly as the wing rotates and the stroke reverses, such that the suction surface of the aerofoil switches sides (figure 5b).

Figure 4

Histograms of summary parameters for the wing kinematics. (a) Wingbeat frequency. (b) Stroke amplitude. The vertical red line shows the kinematics of the standard hovering wingbeat. Horizontal blue lines span the 5th to 95th percentiles for each individual, ranked according to their mean. Note that most of the variation in these parameters is seen between, rather than within, individuals.

Figure 5

Frequency density plots of time-varying wing kinematics over all wingbeats. (a) Wing tip kinematics, showing stroke angle ϕ (red), and deviation angle θ (blue). (b) Blade element kinematics, showing wing pitch angle ω (blue), and aerodynamic angle of attack α (red), both measured at midspan. The shading density corresponds to the frequency density, conditional upon wingbeat phase. Dashed lines indicate ±1 s.d. from the mean; solid lines plot the kinematics of the standard hovering wingbeat.

Aerodynamic force coefficients

With the aerodynamic speed U and angle of attack α defined at the three-quarter chord point, our best estimates for the model parameters were and , such that the lift and drag coefficients (equations (2.14) and (2.15)) may be modelled empirically asandThe variance in the estimated lift and drag curves over 105 random 10% subsamples of the 854 flight sequences was negligible for the lift coefficient but more substantial for the drag coefficient (figure 6b). This reflects the fact that the error in the estimation of the lift coefficient depends only on the error in the estimation of the aerodynamic force coefficient derivative , whereas error in the estimation of the drag coefficient depends also on the error in the estimation of the drag coefficient offset . We therefore investigated the effect of dropping from the regression model in equation (2.23) (i.e. modelling the lift and drag coefficients using equation (2.11) instead of equation (2.15)). This produced a 1% increase in the estimated aerodynamic force coefficient derivative , from 2.70 to 2.73, in compensation for the slight decrease in the predicted drag that would otherwise result from dropping . These changes were associated with only a 0.6% increase in the error sum of squares (table 1), so the inclusion of a drag coefficient offset—though justified theoretically—adds little predictive power to the model. On the other hand, it is clear from figure 6b that the estimated value of the drag coefficient at zero angle of attack is consistently positive across many different subsamples of the data, so the inclusion of in the model is also justified empirically. Also dropping the theoretical added mass correction from the measured forces on the left-hand side of equation (2.23) resulted in an additional 5% increase in the error sum of squares (table 1), which confirms that including the added mass term also has merit in improving the model fit.

Figure 6

Empirically fitted lift and drag polars. (a) Fitted lift coefficient versus fitted drag coefficient across the full range of angles of attack α, for the full model with two empirically fitted parameters (dark blue), and for the reduced model with no drag offset term (cyan). (b) Fitted lift coefficient (orange) and fitted drag coefficient (green) plotted against angle of attack α. The width of the lines indicates their mean ± 1s.d., assessed over 105 random subsamples of the data, each comprising 10% of the recorded flight sequences.

These results assume that the aerodynamic speed U and angle of attack α are defined at the three-quarter chord point. Adjusting the reference point at which the kinematics were defined allowed a small reduction in the error sum of squares, but with the paradoxical result that the estimated drag coefficient offset became negative if the reference point was moved farther than 85% towards the trailing edge (figure 7). Subject to the physical constraint that , the error sum of squares was minimized when U and α were defined at 87% of the chordwise distance back from the leading edge, with and . The chordwise reference point that minimizes the error sum of squares is therefore associated with a binding inequality constraint . This optimization produces only a 1% improvement in the error sum of squares (table 1), and comes at the cost of estimating a third parameter from the data. The best-fitting chordwise reference point cannot be estimated within the existing regression equations, and was instead estimated using an exhaustive search procedure. We therefore prefer to retain the parameter estimates of and with the kinematics defined at the three-quarter chord point in accordance with the prior expectation from classical aerodynamic theory (see §2.5).

Figure 7

Variation in the parameter estimates of the blade element model as a function of the chordwise position of the reference point at which the aerodynamic speed U and angle of attack α are defined. Dark blue line denotes the full model; cyan denotes the model fitted with no drag offset term; note that these models become equivalent beyond the point at which the inequality constraint bites. (a) estimated pressure force coefficient derivative, ; (b) estimated drag coefficient offset, ; (c) associated mean squared error (MSE); see table 1 legend for details. These parameter estimates are made subject to the inequality constraint that , which bites when the kinematics are defined at greater than or equal to 85% chord. See text for further details.

Goodness of fit of the stroke-averaged aerodynamic forces

Because the blade element model is fitted as a regression forced through the origin (equation 2.23), its R2 statistic is not well defined. To assess its goodness of fit with respect to the measured forces, we therefore regressed the fitted aerodynamic forces on the measured aerodynamic forces, without forcing the regression line through the origin (figure 8). We did this separately for each body axis, and found that a large proportion of the variation in the stroke-averaged forces measured in the body’s plane of symmetry was explained in both the x- and the z-axes (R2 = 82.5% and R2 = 79.0%, respectively). By contrast, a much smaller proportion of the measured variation in the stroke-averaged forces was explained in the transverse y-axis (R2 = 18.0%). Moreover, although the regression intercept was appropriately close to zero in all three body axes (x: 0.041 mN; y: 0.000 mN; z: 0.021 mN), the regression slope was only suitably close to one in the body’s plane of symmetry (x: 0.92; z: 1.05), being greatly attenuated by noise in the transverse axis (y: 0.62). This poor fit in the y-axis presumably reflects the fact that the total range of the lateral stroke-averaged aerodynamic forces was small (figure 3e), which leads to a lower signal to noise ratio in the y-axis (figure 8b) than in the x- or z-axis (figure 8a,c). The forces measured in the x- and z-axes are generally well modelled, although a detailed inspection of the regression plots shows that whereas the aerodynamic forces measured in the x-axis are fitted closely over their entire range (figure 8a), the blade element model systematically under-predicts the magnitude of the largest aerodynamic forces produced in the z-axis (bottom left region of figure 8c). Even so, the blade element model does a good job of fitting the overall time-history of the measured stroke-averaged forces on the timescale of an entire flight sequence (figure 9).

Figure 8

Two-dimensional histograms of fitted versus measured stroke-averaged forces for all wingbeats, resolved in the three body axes (a–c). Shading density corresponds to frequency density of data; black line indicates the ideal line of slope one, passing through the origin. Note that at peak aerodynamic force production, the regression model systematically under-predicts the magnitude of the force actually produced in the z-axis.

Figure 9

Time history of the measured and fitted stroke-averaged aerodynamic forces for one flight sequence, resolved in the three body axes (a–c). The forces fitted by the model (red) show a close match to the measured forces (black) in all three body axes. See electronic supplementary material, Video S1 for the wing and body kinematics corresponding to this flight sequence.

Predicted aerodynamic forces through the wingbeat

Because the parameters of the blade-element model were fitted only to the stroke-averaged forces, there is no necessary statistical reason to assume that the resulting model will perform well in predicting the time-varying aerodynamic forces through the wingbeat, but the physical basis of the underlying aerodynamic model is such that it could be expected to. We used the blade element model to predict how the aerodynamic forces are expected to vary through the standard hovering wingbeat that we defined in §2.2 (figures 10 and 11; electronic supplementary material, Videos S2 and S3), to allow us to assess the relative contributions of lift, drag and added mass at different stages of the wingbeat. As a further check on the robustness of our predictions to errors in parameter estimation, we modelled the time-varying aerodynamic forces through this standard hovering wingbeat, across the full range of variation in the aerodynamic force coefficients estimated for the subsamples in §3.3. Despite the variation in the aerodynamic force coefficient parameters fitted in the subsampling analysis (figure 6b), the resulting variation in the predicted time-varying aerodynamic forces was slight in comparison to their variation through the wingbeat (figure 10). The results in figure 10 therefore provide a sound basis for comparing the predictions of our blade element model with future CFD simulations of the standard hovering wingbeat.

Figure 10

Predicted time-varying force components for the standard hovering wingbeat. The plots show the forces predicted for the right wing only, decomposed as lift (orange), drag (green) and added mass (red), together with their sum (blue). (a–c) Forces resolved in the body axes. (d) Vertical aerodynamic force component. The width of the lines indicates their mean ± 1s.d., assessed over 105 random subsamples of the data, each comprising 10% of the recorded flight sequences.

Figure 11

Spanwise distribution of predicted aerodynamic force at four phases of the standard hovering wingbeat. (a–d) Resultant aerodynamic force (blue arrows). (e–h) Quasi-steady lift (orange arrows) and drag (green arrows) force, together with unsteady added mass force (red). The four phases of the wingbeat correspond to (a,e) beginning of the downstroke; (b,f) middle of the downstroke; (c,g) beginning of the upstroke; (d,h) middle of the upstroke. The hoverfly body is shown in the background, and the grey line shows the wing tip path. The viewpoint rotates with the stroke, so as to look down the spanwise axis of the wing at all times. Individual blade elements are plotted as black lines, to show the spanwise distribution of wing twist through the wingbeat. See electronic supplementary material, Videos S2 and S3 for animated versions of these forces at all time steps.

Because the stroke plane is close to horizontal during the standard hovering wingbeat, weight support is attributable primarily to aerodynamic lift, which is predicted to account for 77.2% of the stroke-averaged vertical force (figure 10d). The added mass and drag forces each make a small net positive contribution to weight support, providing 14.3% and 8.5%, respectively, of the predicted stroke-averaged vertical force. The predicted lift force peaks at the middle of each half-stroke, when the wing’s translational velocity is highest (figure 11f,h; electronic supplementary material, Video S3), so this is also the phase of the stroke during which the majority of the vertical force is expected to be produced (figures 10d and 11b,d). By contrast, the added mass force peaks at the beginning of each half-stroke, when the wing’s acceleration normal to its chord is highest (figure 11e,g; electronic supplementary material, Video S3). The drag force has a more complicated time-history again, peaking at several points through the stroke (figure 11e–h; electronic supplementary material, Video S3). Aerodynamic force production in the transverse y-axis is qualitatively similar on both the upstroke and downstroke (figure 10b), and the same is true of the vertical component of the aerodynamic force (figure 10d), but in each case the amplitude of the forces is somewhat diminished on the upstroke relative to the downstroke. By contrast, aerodynamic force production along the x- and z-axes of the body displays a marked asymmetry between the upstroke and the downstroke (figure 10a,c). The dynamics of hovering force production are therefore considerably more complex than might first appear from the time-history of the vertical aerodynamic force component alone (figure 10d).

Discussion

Perhaps the greatest challenge in modelling insect flight is to predict the aerodynamic forces that the flapping wings impart as they undergo a variety of complex aeroelastic motions in a variety of different flight conditions. Although current CFD techniques allow these forces to be predicted with great accuracy for a given set of wing or body kinematics, it remains extremely time consuming to compute the flows associated with variable wing or body kinematics or with different wing morphologies [5-8]. Moreover, the aerodynamic assumptions that classical analytical models must make to fix a solution for the aerodynamic force coefficients are so restrictive as to prevent their realistic application to insect flight [12]. Here we have aimed to find a middle ground, by estimating the aerodynamic force coefficients empirically for a simple analytical blade element model that captures the scaling of the forces expected from first principles, and which models the forces with sufficient accuracy to be used in a range of other applications. By estimating just two numerical parameters—the derivative of the pressure force coefficient () and the drag coefficient offset ()—for the largest dataset of insect wing and body kinematics obtained to date [18], we have captured 80% of the total variation in the measured stroke-averaged aerodynamic forces in the sagittal plane over N = 26 541 wingbeat pairs recorded in freely hovering and manoeuvring Eristalis hoverflies. This is possible because whereas the aerodynamic form of the model is quite simple and uses empirical estimates of the force coefficients that are averaged over many different flight conditions, the modelling of the kinematics is sufficient to capture almost the full complexity of the wingbeat. In summary, the key contribution of this paper is to provide a kinematically accurate blade-element model of the stroke-averaged forces of insect flight, fitted and validated with respect to an extensive free-flight dataset including a wide range of flight manoeuvres.

Key features of the modelling

A key feature of our kinematic model is that it takes full account of the torsional deformation of the wing, without either limiting the motion to a fixed stroke plane, or treating the entire wing as a flat plate. This matters, because wing flexion is a defining characteristic of insect flight, which can reduce its aerodynamic power requirements and enhance the useful force produced [8,44]. For example, when we tried treating the wing as a flat plate operating at a pitch angle matched to that of the twisted wing at mid-span, we found that the error sum of squares increased by 30%, which was associated with a large decrease in the estimated value of and an unrealistically large increase in the estimated value of (table 1; electronic supplementary material, figure S3c). Modelling the kinematics accurately is therefore at least as important as modelling the flow accurately, so it does not make sense, for example, to go the effort of solving the full Navier–Stokes equations numerically for an oversimplified model of an insect’s wing kinematics.

Another key feature of the kinematic modelling is that it implicitly accounts for the body’s own velocity and acceleration when computing the velocity, angle of attack and acceleration of each blade element. The body’s motion plays a key role in flight stability (e.g. [45,46]), so it is essential that any aerodynamic model which aims to investigate free-flight behaviour takes account of these effects. Moreover, whereas the velocity of the body is usually considerably smaller than the velocity of the wing tip, the two can become comparable in magnitude during fast manoeuvres and fast forward flight. In fact, there is clear evidence of the importance of the body’s motion in our modelling, because if we ignored the body’s velocity and acceleration when fitting the blade element model, then the error sum of squares increased by 34%, and the estimate of the drag coefficient offset became unreasonably high (table 1; electronic supplementary material, figure S3d).

Concerning the aerodynamic modelling, it is noteworthy that moving the reference point at which the kinematics are defined backwards from its default three-quarter chord position results in at best a 1% reduction in the error sum of squares and is associated with the drag coefficient offset being driven unrealistically towards zero. Conversely, moving the reference point forwards from the three-quarter chord point causes a rapid worsening of the model fit, which implies that the default approach is successfully capturing the rotational lift as expected (see §2.5). This provides a useful empirical validation of the approach of treating rotational lift together with translational lift [16] by defining the aerodynamic speed U and angle of attack α at the three-quarter chord point [12], consistent with the results of previous work using mechanical flappers [7].

Comparison to other models

Blade element models have been widely used to predict the aerodynamic forces of insect flight, but never before in combination with such detailed kinematic data from free-flying insects. Our model differs from others used previously in the following respects: (i) it includes the full three-dimensional motion of the wing, including torsional deformation; (ii) it incorporates all six rotational and translational degrees of freedom of the body; (iii) it models the rotational aerodynamic force by calculating the angle of attack at the three-quarter chord point rather than treating this as a separate contribution; (iv) it captures systematic variation in the direction of the quasi-steady aerodynamic force with respect to the wing’s surface and the flow; and (v) it has a simple and transparent aerodynamic form developed from physical principles. The first two of these distinguishing features relate partly to the availability of data, but all are fundamental.

Concerning our treatment of the rotational forces, most recent blade element models include a separate rotational lift component, sometimes called a rotational circulation force (e.g. [7,15]). Indeed, the blade element model of Nakata et al. [14], with coefficients fitted using computational fluid dynamics also includes a separate rotational drag force. These rotational forces are expressly intended to capture the aerodynamic effects of pitching or twisting about the spanwise axis of a flapping wing, but by defining the velocity of each blade element at its three-quarter chord point [12], all or most of these rotational effects can be incorporated directly into the calculations of translational lift and drag, thereby simplifying the model conceptually and reducing the number of free parameters that must be estimated [16].

The direction of the quasi-steady pressure force is determined in our model by the balance of the orthogonal lift and drag forces, where the lift on each blade element is assumed to act perpendicular to the relative air flow, and where the drag is assumed to act in the direction of the relative air flow. Other blade element models have defined the circulatory force as acting normal to the blade element, but have assumed that its magnitude is equal to that of the resultant lift and drag [6,7,9,15], which means that a circulatory force is assumed to act perpendicular to the chord even at zero angle of attack if there is any friction drag or pressure drag under these conditions. By contrast, the treatment of the drag offset term in our model means that the direction of the quasi-steady aerodynamic force varies with respect to the chord (figure 12a). Specifically, whereas the quasi-steady force is almost perpendicular to the chord at angles of attack greater than about 45°, it becomes tilted back substantially at lower angles of attack, ultimately becoming tangent to the chord at zero angle of attack.

Figure 12

Predicted direction and relative contribution of the quasi-steady aerodynamic force across different angles of attack. (a) The direction of the resultant quasi-steady force vector is represented by its angle (γ) with respect to the blade element chord, across the full range of measured aerodynamic angle of attack (α). The predicted angle γ is calculated as the angle between the resultant lift and drag force and the blade element (inset diagram), using the regression estimates for and . This angle γ only approaches 90° (dotted line) at angles of attack above approximately 45°. (b) Distribution of predicted contribution to the stroke-averaged quasi-steady aerodynamic force at different aerodynamic angles of attack for all recorded wingbeats combined. The horizontal red bars display the total percentage of the force produced at angles of attack above and below 45° (i.e. the area under the graph to either side of the vertical red line).

This behaviour is appropriate given the inevitable presence of friction drag, and the likely importance of pressure drag even at low angles of attack. Although insect wings are commonly treated as approximating an idealized thin aerofoil, they do in fact have a finite thickness, particularly at the leading edge, which typically functions as a reinforced spar. Wing corrugation due to venation also makes wings inherently three-dimensional structures, which can further increase profile drag [10,11]. Aerodynamic measurements of real insect wings [47,48], or mechanical and computational models thereof [3,49], have all indicated the presence of significant drag at zero angle of attack, and hence of the tilting back of the resultant aerodynamic force at low angles of attack. While it is true that insect wings operate at characteristically high angles of attack, in the hoverfly data that we have presented here, 47% of all of the aerodynamic force produced by the wing is predicted to occur at angles of attack α < 45° (figure 12b), with a mode at α = 36°, for which the force vector will be tilted back significantly from 90° (figure 12a).

Perhaps the most important feature of our model is the simplicity of its form, which despite being highly nonlinear in the kinematics, is nevertheless linear in its two free parameters and . This means that we are able to fit these parameters analytically using linear least squares, guaranteeing that they are optimized globally and at speed—even on the very large dataset that we employ. Moreover, because the form of our model was developed from first principles rather than with reference to the data, it provides useful physical insight into the ad hoc form of some influential models that have been fitted previously. In particular, Dickinson et al. [3] modelled the lift and drag coefficients for their robotic flapper as C = 0.225 + 1.58 sin(2.13α − 0.14) and C = 1.92 − 1.55 cos(2.04α − 0.17), where α is in radians and where all of the numerical constants are estimated from the data. These complicated formulae use a total of eight free numerical parameters to model C and C, but make better sense physically when it is observed that they approximate trigonometric double angle formulae. Noting that sin(2α) = 2 sin α cos α and cos(2α) = 1 − 2 sin2α, equations (3.1) and (3.2) describing our own fitted model can be restated asandwhich are not dissimilar to the equations fitted by Dickinson et al. [3] that themselves approximate trigonometric double angle formulae. This in turn reflects the fact that whereas the resultant force varies sinusoidally with the angle of attack, its direction is approximately normal to the wing’s surface, such that its lift and drag components resolved perpendicular and tangent to the relative airflow vary as cosine and sine functions of the angle of attack, respectively.

Extension to other datasets

The same blade element model can be straightforwardly applied to other insect species, because other than the morphological parameters of wing length and wing shape (figure 2), there are no modelling assumptions that are specific either to hoverflies or to the dataset that we used. Indeed, one of the strengths of our approach is that the same blade element model can be straightforwardly fitted to any similar dataset, using linear least squares to optimize the force coefficient parameters: Matlab code implementing the blade element model is provided as Supporting Data S1 to this end. It is clear from our modelling that knowledge of the wing twist distribution is essential to fitting the forces accurately (table 1), which reinforces the need for future kinematic studies to measure wing deformation as we have done here. Of course, in design problems where the wing twist distribution is a parameter that may need to be optimized rather than measured, the inclusion of wing twist in our analytical blade element model makes it suitable for use in fast global optimization of the wing deformation parameters prior to local refinement using CFD.

Our analysis shows that the blade element model is robust to drastic reductions in sample size, with data subsampling producing comparatively small changes in the predicted forces when fitting the force coefficients to random subsamples comprising only 10% of the recorded data. This means that the same modelling approach can be applied to datasets much smaller than the N = 26 541 wingbeats that we analyse here. If necessary, the simplicity of the model can be increased further by excluding the small drag offset term. This then means that only the derivative of the pressure force coefficient need be estimated, which reduces the variance of the parameter estimate (electronic supplementary material, figure S3b). Given the comparatively weak signature of the drag offset term at high angles of attack, this may be a preferable approach when working with smaller or noisier datasets than the one available here. However, while the inclusion of the drag offset term provides only a 0.6% reduction in the error sum of squares, we retain it in our model because its importance is supported by both theory and experiment [3,49]. Furthermore, although its inclusion has a minimal effect on the accuracy of our modelling of this dataset for Eristalis, the drag offset term may be more important for other species with different morphologies and especially for those operating at lower Reynolds numbers.

As with any form of regression modelling, an important caveat is that the dataset must contain sufficient variation in the independent variables to enable a good fit (compare, figure 8a,c with figure 8b). In principle, our regression estimates of the lift and drag coefficients could be replaced by estimates from model wings [3] or CFD [14]. However, an obvious risk of this approach is that the aerodynamic properties of a model wing, or even those of a detached wing suffering rapid desiccation [50], may differ markedly from the aerodynamic properties of a real wing in vivo. These problems are avoided completely by our approach of fitting the parameters of the blade element model empirically to free-flight data from live insects.

Limitations

Although our blade element model fits the aerodynamic forces well in the sagittal plane over most of their range, it systematically under-predicts the magnitude of the largest aerodynamic forces produced in the z-axis (figure 8c). This discrepancy presumably indicates a nonlinearity in aerodynamic force production that the quasi-steady blade element model fails to capture. More specifically, as most of the force in the z-axis is produced on the downstroke (figure 10c), we hypothesize that this nonlinearity reflects some unsteady aerodynamic mechanism that becomes increasingly important as force production is increased on the downstroke. One obvious possibility is that this nonlinearity is due to the presence of a leading-edge vortex (LEV) on the functional upper surface of the wing. This is one of the characteristic aerodynamic mechanisms of insect flight, allowing the lift curve to be extended to high angles of attack by delaying stall [51]. Delayed stall is already implicit in our model because of its explicit assumption that the wing does not suffer a sudden loss of lift at high angles of attack (figure 6). It remains unclear, however, whether the presence of an LEV enhances the aerodynamic force coefficients by amplifying some portion of the lift curve rather than merely by extending it [52]. Hence, although the average aerodynamic effect of the LEV should be captured by the empirical force coefficients that we have estimated from our data, it is plausible that the model might still underestimate the lift enhancement provided by the LEV at very high angles of attack. It is also worth noting that our model does not account explicitly for the three-dimensional effects of spanwise flow impacting the strength and stability of the LEV [13], nor for the possible effects of interactions between the wing and the wake shed on the preceding half-stroke [7].

A second possibility is the clap-and-fling mechanism [12], which can occur if the wings approach one another closely at the top of the upstroke and are then flung apart on the downstroke (figure 13a). This mechanism may also be modified by the effects of spanwise bending and chordwise camber [53], which we do not model directly here. Interestingly, although there is a positive association between the stroke amplitude and the magnitude of the measured aerodynamic force along both the x- and z-axes (figure 13b), as expected under a quasi-steady model of the forces, this association is much stronger in x (R2 = 0.53) than in z (R2 = 0.27). This suggests that increases in the magnitude of the aerodynamic forces in z are not principally driven by increases in stroke amplitude. On the other hand, there is a negative association between the wing tip separation at the start of the downstroke, and the magnitude of the measured aerodynamic forces (figure 13c), which is stronger in the z-axis (R2 = 0.62) than in the x-axis (R2 = 0.25). This negative association would be expected under an unsteady clap-and-fling mechanism, and the strength of the association in z suggests that increases in the magnitude of the aerodynamic forces in this axis may indeed by driven by an unsteady clap-and-fling mechanism. This is consistent with the interpretation that wing-wing interactions affect force production in the z-axis more than in the x-axis. Since the blade element model will not capture this nonlinearity directly, it follows that the linear parameter estimate for may overestimate the true value of the quasi-steady force coefficient derivative . Another possible reason for this discrepancy is in the modelling of the added mass forces, where there remains room for further improvement as detailed in the Methods.

Figure 13

Histograms of wing tip separation and stroke amplitude, and their association with the measured stroke-averaged aerodynamic force. (a) Histograms of wing tip separation at the start of the downstroke (shaded bars) and at the start of the upstroke (unshaded bars). (b,c), Frequency density plots showing stroke amplitude and wing tip separation at the start of the downstroke, versus the measured forces in the x (blue) and z (red) axes of the body. Shading density corresponds to frequency density of data. Wing tip separation was calculated as the distance between the wing tips at the start of the half-stroke, normalized by the mean wing chord; R2 statistics for the linear associations between the wingbeat parameters and the measured forces are shown for each axis.

Finally, although the effects of body motion are captured in our modelling of the wing kinematics and aerodynamics, it is important to note the body itself will produce drag—and perhaps some lift—in forward flight [48]. Modelling the aerodynamic forces produced by a bluff body is not straightforward, owing to the likelihood of sudden flow separation above some critical angle of attack, but as most of the flight sequences that we modelled were close to hover, it is reasonable to assume that the aerodynamic forces on the body would have been overwhelmed by the aerodynamic forces acting on the wings at most stages of the wingbeat.

Conclusion

We have shown here that an analytical blade element model with just two empirically fitted coefficients provides a close fit to the measured stroke-averaged aerodynamic forces of free-flying insects in manoeuvring flight. The alternative approach of using computational fluid dynamics modelling is capable of capturing fine aerodynamic detail, and has even led to the discovery recently of novel unsteady mechanisms (e.g. [1]), but is computationally expensive, taking many orders of magnitude longer to deliver results than the analytical blade element model presented here. Both approaches therefore have a complementary role to play. Analytical blade element modelling functions well for investigating large datasets, studying the effect of changing wing kinematics parametrically, and making quick comparisons across species. Conversely, a numerical approach is preferable where high-fidelity predictions, detailed time histories, or insight into unsteady aerodynamic mechanisms is required. The strengths of the analytical blade element model that we have presented here are (i) the simplicity of its underlying aerodynamic equations; (ii) the complexity of the deforming wing kinematics and body motions that it models; and (iii) the fact that its aerodynamic force coefficients are fitted empirically to free-flight data from real insects, thereby capturing the full scope of the insect’s flight dynamics. Besides demonstrating the importance of prioritizing accurate modelling of the deforming wing kinematics ahead of detailed modelling of the fluid dynamics, we expect that our model will serve as a useful, validated tool for future research on insect flight dynamics and control.