A semi-empirical aerodynamic model of insect flight forces

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

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 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 [4]. Dipteran flies have received particular attention, because their possession of only two functional wings halves 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 utilised, each with their own advantages and disadvantages. More sophisticated techniques include the use of mechanical flappers [5][6][7] or computational fluid dynamics [e.g. 8,15,24,33]. Both approaches determine the aerodynamic forces from a predefined set of wing kinematics, allowing the effect of specific kinematic parameters to be investigated experimentally [7,23,33], 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 [28,30]. 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 [11]. To the extent that the flow around a real wing is inherently three-dimensional and coupled to the wake, 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 summarise 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 [21,24,26], 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 [24]. 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 [21,24,26]. The key is to identify an appropriate analytical formulation describing how the aerodynamic forces on each blade element vary with the wing kinematics, 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 (15,21,24,26; but see 27). In each case, these models have been parameterised by a set of empirical force coefficients fitted to measurements made under different kinematic conditions using either a mechanical flapper [15,24,26] or a numerical model [21]. The resulting blade element models involve either one [15,24,26] or two [21] 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 [21] to four [24,26] or five [15] 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 [21] to nine [24,26] or even eleven [15] empirically-fitted parameters. This over-parameterization 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 [15,21,24,26], it is unknown how well they predict the aerodynamic forces and moments on real insects undergoing free-flight manoeuvres involving complex wing deformations.
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 empirically-fitted free parameters. We achieve this by using linear least squares modelling to fit the numerical coefficients of an unsteady blade element model 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 make sense of some of the more complicated empirical functions that have been fitted to model the aerodynamic force coefficients previously [6,15] and that have been co-opted into subsequent models [24,26]. 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. Whilst 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 [32], 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. Indeed, for the largest insects, in which the body dynamics operate on a similar timescale to the wingbeat [32], it should even be possible to fit the aerodynamic force coefficients directly to the time-varying rather than stroke-averaged forces.

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 [31]. For completeness, we begin by providing a brief summary of the experimental methods and kinematic reconstruction technique, before providing a detailed description of the aerodynamic modelling that forms the primary contribution of this paper.

(a) 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 synchronised high-speed video cameras (SA3, Photron Ltd, Bucks, UK) with 180mm 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 synchronised 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 was far beyond the range of the visible spectrum for Eristalis [upper limit: 600 nm; 1,2,16], 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 2D calibration grid held in a range of positions and orientations [29]. We captured between 10 and 50 separate recordings of each hoverfly, before anaesthetising it with CO 2 at the end of the experiment, and weighing it using a microbalance with 0.1µg readability (UMX2, Mettler Toledo Ltd, Leicester, UK). We analysed all of the flight sequences in which both wings were visible to all four cameras for ≥ 5 wingbeats, 4 giving a total sample of 854 flight sequences from n = 36 hoverflies. A fully automated shapecarving procedure was used to label voxels contained within the minimum convex hulls of the insect's body and its two wing outlines, respectively [31].

(b) Flight kinematics modelling
Throughout the paper, we use boldface symbols to represent vector quantities, and use t to represent time. We defined the body kinematics using a right-handed, rotating, body-fixed axis system {x b , y b , z b } with its origin at the centre of volume of the body voxels, its x b -axis directed anteriorly along the major axis of the body voxels, and its y b -axis pointing rightward, parallel to the line connecting the wing bases (Fig. 1). We measured the position of the body axes, X b (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 smoothing spline with a tolerance chosen to give the same degree of smoothing as a 3rd order Butterworth filter with a -3 dB cut-off frequency of 100 Hz (i.e well below the insect's 188 ± 14 Hz wingbeat frequency; mean ± SD). We doubledifferentiated each spline analytically to estimate the insect's instantaneous linear acceleration with respect to the laboratory coordinate system,Ẍ b (t), resolved in the lab axes. with its origin at the centre of volume of the insect's body (grey silhouette), its x b -axis directed anteriorly along the major axis of the body voxels, and its y b -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 R , y R , z R } (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 L , y L , z L } (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 R -axis, then through its deviation angle −θ about the x R -axis, so as to bring the line from wing base to wing tip into alignment with the y R -axis. The pitch angle ω was then defined as the angle from the x R y R -plane to the anatomical ventral side of the rotated chord perpendicular to the y R -axis. (c) The speed of a blade element (U ) is measured at its three-quarter chord point (green dot), and is defined as the hypotenuse of the components of the velocity vector directed parallel to the blade element chord (U c ) and normal to the blade element surface (U ⊥S ). 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 R , y R , z R } parallel to the rotating body axes {x b , y b , z b }, 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 5 wing ( Fig. 1). We defined the local pitch angle ω(r) of the wing at radial coordinate r by rotating the wing's measured outline through the angle −φ about z R , then through the angle −θ about x R , so as to bring its spanwise axis into alignment with y R . The local pitch angle ω(r) of the wing was then defined as the angle from the x R y R -plane to the anatomical ventral side of the chord, perpendicular to the y R -axis at radial coordinate r (Fig. 1B). We summarized the spanwise twist of the wing by regressing ω(r) on r, modelling the local pitch angle as: where ω 0 is the pitch angle offset, and ωr is the linear twist gradient. The kinematics of the left wing were defined independently using an equivalent left-handed axis system {x L , y L , z L }.
We filtered our measurements of θ, φ, ωr, and ω 0 for each wing using the same quintic spline method as we used to smooth the body kinematics, matching the tolerance factor to a cutoff frequency of 500 Hz for the tip kinematics θ and φ, and 800 Hz for the twist kinematics ωr and ω 0 . The kinematic data were 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 wingbeats of differing period, we used cubic interpolation to resample the wing and body kinematic data and total aerodynamic force at 100 evenly-spaced time steps through each wingbeat, beginning at the start of the downstroke. On average, this means that the kinematic measurements were up-sampled by a factor of approximately 5.

(c) Flight dynamics modelling
The goal of this paper is to use aerodynamic modelling to relate free-flight measurements of body kinematics to the wing kinematics that produce them. Beginning from first principles, the only ways in which a fluid can impart force to a solid surface are through pressure forces acting normal to the surface, and through 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 as: where m is the insect's mass,Ẍ b (t) is the acceleration of the insect's body with respect to the laboratory coordinate system, g is gravitational acceleration, P is the total pressure force, and F is the total friction force. All of these vector quantities are resolved in the non-rotating axes defined by the lab-fixed coordinate system. Since the wingbeat period of a hoverfly is much shorter than any characteristic timescale of its body's dynamics [32], 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: where n is wingbeat number. This 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 [32], it might be more appropriate to retain the original form of Eq. 2.2. The lefthand side of Eq. 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 estimate of the time-varying aerodynamic force that we make later using our blade element model. The stroke-averaged expression in Eq. 2.3 forms the basis of the empirical estimation of the aerodynamic forces in this paper.

(d) Aerodynamic modelling
We begin this section by developing the scaling of the aerodynamic forces, to motivate the theoretical form of the kinematic predictor variables that we use to fit the aerodynamic coefficients of the quasi-steady blade element model, and to introduce the theoretical form of the equations that we use to model the unsteady forces.
(i) Functional form of the quasi-steady aerodynamic forces The ratio of inertial to frictional forces on a flapping wing is given by its chord Reynolds number, which is of order 10 3 for Eristalis. This implies that the pressure forces will dominate the friction forces in Eqs. 2.2 and 2.3, which in turn reflects the fact that the shear layer surrounding the wing is thin in relation to its chord. It follows that viscous shear in this inner region of the flow cannot ultimately be responsible for setting fluid into motion in the outer region of the flow. The flow in this outer region must instead be driven by the pressure gradients that result ultimately from the imposition of the boundary condition that fluid cannot cross the solid surface of the wing. This boundary condition implies that at every point on the wing's surface, the surface-normal component of the flow (Q ⊥S ) must be the same as the surface-normal component of the wing's own velocity (V ⊥S ), which we may write locally as: where β is the local incidence angle between the surface of the wing and its local velocity, and where V is the local speed of the wing. The pressure gradients resulting from the global imposition of this boundary condition will in general modify the tangential flow at the wing's surface, resulting in a net circulation (Γ ) defined as the line-integral of the tangential flow around the airfoil. The strength of the resulting circulation cannot be determined without making further and more restrictive assumptions, but must obviously depend on the magnitude of the normal flow in Eq. 2.4. For a flat plate airfoil undergoing pure translational motion through still air, the local incidence angle β is the same everywhere on the upper and lower surfaces of the wing, and is equal to the overall angle of attack α between the velocity of the airfoil and its chord. If the airfoil is also rotating steadily, then the instantaneous effect of any nose-down rotation, say, will be to increase β ahead of the axis of rotation, and to decrease β behind it. We must therefore specify how to weight these local contributions to the overall angle of attack α, which for a linear weighting distribution is equivalent to specifying a unique point on the airfoil at which the boundary condition in Eq. 2.4 is to be satisfied. Under classical inviscid lifting-line theory, this collocation point falls three-quarters of the way back from the leading-edge if the flow is assumed to depart smoothly at the trailing edge. The same collocation point also appears in the inviscid theory of the smallamplitude unsteady motion of a pitching and plunging flat plate airfoil, so for a flat plate airfoil undergoing quite general motion in a potential flow, the quasi-steady effects of wing rotation can be adequately captured by treating the overall angle of attack α as being equal to the local incidence angle β at the three-quarter chord point [11].
The form of the boundary condition in Eq. 2.4 therefore suggests that under steady or quasisteady conditions, the circulation Γ around an airfoil can be expected to vary as a function of cU sin α, where c is the chord length, and where the overall speed U and angle of attack α of the airfoil are both assumed to be measured at its three-quarter chord point. Noting that the Kutta-Joukowski theorem for two-dimensional inviscid flows predicts that the pressure force varies as ρU Γ where ρ is the fluid density, we see that this inference is consistent with the classical result that the quasi-steady pressure force is in fact fixed at a value πρU 2 c sin α per unit span for a flat-plate airfoil operating at low enough angle of attack that the flow departs smoothly from the trailing edge [17]. Whilst there is no necessary expectation of linearity in sin α under other flow conditions, any aerodynamic model that we fit must naturally be able to handle this classical 7 potential flow regime as well as the more realistic separated flow conditions applicable to insect flight. This being so, it is reasonable to attempt to model the quasi-steady pressure force (Pqs) on a blade element of width ∆r using the scaling: where the unknown constant of proportionality and direction of action of the force remain to be determined empirically.
The two-dimensional Kutta-Joukowski theorem predicts that this quasi-steady pressure force acts perpendicular to the relative velocity of the airfoil at low angles of attack, but some departure from this is inevitable with the separated flows experienced at higher angles of attack on a threedimensional wing. In fact, the pressure force acts more nearly perpendicular to the chord when the flow is separated over the upper surface of the airfoil [17]. This contradiction is not too problematic in practice, because even at low angles of attack, the flow induced by the wake has the effect of tilting the aerodynamic force vector back considerably on a low aspect ratio wing, bringing its line of action more nearly perpendicular to the chord. In any case, given that the stroke-averaged pressure force is inevitably dominated by the very high pressure forces experienced during phases of the wingbeat when the wing is at very high angle of attack, it is reasonable to assume-as a first approximation-that the quasi-steady pressure force acts normal to the chord, signed positive in the direction of increasing angle of attack. Decomposing this force into a quasi-steady lift component, L = Pqs cos α, and a quasi-steady drag component, D = Pqs sin α, then yields the scaling: where the constant of proportionality is assumed to be the same in both expressions. Note that, by definition, the lift acts perpendicular to the relative flow, whereas the drag acts in the direction of the relative flow.
(ii) Quasi-steady force coefficients By convention, the non-dimensional lift coefficient (C L ) and drag coefficient (C D ) of a blade element are defined by the identities: and similarly for the pressure force coefficient (C P ) with respect to the overall quasi-steady pressure force Pqs. Combining these identities with the scalings in Eqs. 2.6-2.7 implies that the local lift and drag coefficients may be approximated as: where the unknown constant C Pα remains to be determined empirically, but can be interpreted 8 as the derivative of the pressure force coefficient with respect to the angle of attack, at vanishing angle of attack.
Although the pressure forces are expected to dominate the friction forces at Reynolds numbers of order 10 3 , this does not mean that the friction forces can be dismissed entirely in the drag force calculations-especially at low angles of attack, for which the pressure forces are expected to be small. Eq. 2.11 is missing the effects of friction drag, but as friction drag depends only on the tangential flow in the boundary layer flow, it does not depend strongly on the angle of attack. The simplest method of accounting for it is therefore by adding a drag offset term to the righthand side of Eq. 2.11: Here, the unknown parameter C D0 is simply the drag coefficient at zero angle of attack, which includes the friction drag, together with any non-zero pressure drag that may be present at zero angle of attack, and which again remains to be determined empirically.

(iii) Functional form of the unsteady forces
The quasi-steady aerodynamic force must be supplemented by an unsteady aerodynamic force whenever the wing accelerates, as can be shown by rewriting the normal component of the local flow velocity in Eq. 2.4 as V ⊥S = V sin β, and differentiating this restated boundary condition with respect to time t to yield: which shows that any acceleration of the wing in the direction normal to its surface must be accompanied by an equal acceleration of the fluid at the wing's surface. In accordance with Newton's laws of motion, the reaction to this acceleration of the fluid must be an unsteady pressure force (Pus) acting in the opposite direction normal to the wing's surface. Naturally, the surrounding fluid must also be accelerated, to a degree that will decline away from the wing, but the strength of the reaction force can be expected to scale linearly with the normal component of the wing's acceleration, so that the resulting forces will vary as if the wing were accelerating a virtual mass of fluid with it. For this reason, the unsteady pressure force is usually referred to as an added mass force. The added mass force is complicated to model for a body of arbitrary shape undergoing arbitrary motion, but throughout the insect flight literature, it is conventionally approximated using the classical potential flow solution for a flat plate in translational motion at a constant rate of acceleration. This gives the the unsteady pressure force for any given blade element as: whereU ⊥S is the acceleration of the blade element normal to its surface, defined as the projection of its acceleration onto the surface normal, signed positive in the direction of increasing angle of attack. Note that whereas the velocity of a blade element is defined at the three-quarter chord point for the purposes of determining the quasi-steady aerodynamic forces, the acceleration of the blade element is defined at the half-chord point for the purposes of determining the unsteady added mass forces. Unlike the inertial forces associated with the wing's own acceleration, the added mass force will not sum to zero over a periodic flapping cycle, except under some specific and generally unrealistic symmetry conditions. Nevertheless, the mean of the added mass force will be small in comparison to the scale of its variation through the wingbeat, and on grounds of low signal to noise ratio, it therefore makes sense to model it directly using Eq. 2.14, rather than to attempt to model it empirically using the regression technique that we use for the quasi-steady aerodynamic forces in the next section.
(iv) Blade element modelling The final step in the aerodynamic modelling is to use the equations in Section ii 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 be used together with Eqs. 2.2 and 2.14 to formulate a set of linear equations that can be solved for the unknown aerodynamic force coefficients C Pα and C D0 , 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 (Fig. 2). We may therefore model the instantaneous lift (L), drag (D), and added mass (A) forces acting on the ith blade element at time t as: in which 1 U , 1 ⊥U,r , 1 ⊥S are unit vectors resolved in the body axes and defined as follows. The unit vector 1 ⊥U,r 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 1 U 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 ⊥S 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 1 ⊥S = 1 ⊥U,r cos α + 1 U sin α.
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 Supplementary Methods, such that Eqs. 2.15 and 2.16 account implicitly for the effect of wing rotation on the quasi-steady pressure force (see Section i). In contrast, the normal accelerationU ⊥S is defined at the half-chord point for the purposes of determining the added mass force (see Supplementary Methods). All of these quantities measure the blade element kinematics with respect to the laboratory coordinate system, and therefore account implicitly for the effects of the body's motion during forward flight and manoeuvring. Our kinematic measurements record the motion of each wing at 100 discrete time steps 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 1 ⊥S = 1 ⊥U,r cos α + 1 U sin α to eliminate two of the trigonometric terms when combining Eqs. 2.15 and 2.16. By an obvious use of notation for the summation terms, we will abbreviate the righthand side of Eq. 2.18 as C Pα Σ Pα + C D0 Σ D0 + Σ A . It follows that the total stroke-averaged aerodynamic force on both wings may be modelled for each wingbeat as: where the superscripts L and R refer to the summations for the left and right wings, respectively, and in which it is assumed that these summations are resolved in the lab axes for consistency with the definition of the terms on the lefthand side. Substituting this expression in to Eq. 2.3 for the measured stroke-averaged aerodynamic force and rearranging, we have the expression: for each wingbeat, whereˆ = [ˆ X ,ˆ Y ,ˆ Z ] is an error term that accounts for the approximation in the modelling in each of the three lab axes. Eq. 2.20 is linear in the unknown force coefficient parameters C Pα and C D0 , and can therefore be solved using linear regression to provide parameter estimatesĈ Pα andĈ D0 that minimise the sum of the squared error over all n ∈ 1 . . . N wingbeats and over all three axes. This approach allows us to estimate the quasi-steady lift and drag coefficients of hoverflies empirically by taking full advantage of the very large sample size of N = 26, 541 wingbeats available in our dataset.

Results (a) 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 (Supplementary Video S1) includes brief periods of slow forward flight, punctuated by 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 [12,13,18]. This range of behaviour is reflected in the variability of the measured aerodynamic forces acting along the x b -axis (0.97 ± 0.43 mN) and z b -axis (-1.11 ± 0.4 mN) of the body (mean ± SD; Fig. 3a,c). The forces measured along the y b -axis were much less variable (0.00 ± 0.11 mN; Fig. 3b), consistent with the orthodoxy that comparatively little lateral aerodynamic force is produced during manoeuvres [3]. Nevertheless, the double derivatives of body position vary to a similar extent in all three body axes (Fig. 3d-f), indicating that the asymmetry of force production in the x b -and z b -axes is explained by the need to overcome the body's acceleration due to gravity, which is small in the y b -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. (d-f ) Measured stroke-averaged acceleration,Ẍ b (n), 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 b -and z b -axes of the body.

(b) Wing kinematics
Wingbeat frequency and stroke amplitude both vary greatly over the dataset, but their variability owes more to variation between individuals than within (Fig. 4), and the time-history of the wing kinematics is actually quite stereotyped over the whole dataset (Fig. 5B). Because the wing tip trajectory is inclined at approximately 45 deg to the body, the stroke angle φ and deviation angle θ always have similar oscillation amplitudes, and both vary approximately sinusoidally through the wingbeat (Fig. 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 aerodynamic suction surface of the airfoil switches sides (Fig. 5B). In light of the very large sample size, and because the regression model does not take account of autocorrelation in the forces from one wingbeat to the next, we do not report 95% confidence intervals for these parameter estimates. Instead, we test the robustness of the analysis using subsampling, repeating the regression modelling 10 5 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 estimates of the aerodynamic force coefficients that results from the variation between different individuals and different flight sequences, and does so at a sample size that is more realistic for future studies than the exceptionally large sample used here (approximately 10 2 rather than 10 3 flight sequences). The variance in the lift and drag curves estimated with random 10% subsampling is negligible in respect of the lift coefficient, but is more substantial in respect of the drag coefficient ( Fig  6). This reflects the fact that error in the estimation of the lift coefficient depends only on the error in the estimation of the aerodynamic force coefficient derivative C Pα , whereas error in the estimation of the drag coefficient additionally depends on the error in the estimation of the drag coefficient offset C D0 . We therefore investigated the effect of dropping the drag offset term C D0 from the model (Eq. 2.20). If C D0 is set to zero, then the estimated aerodynamic force coefficient derivativeĈ Pα increases slightly, from 2.70 to 2.73, thereby compensating for the drop in predicted drag that would otherwise result from removing C D0 . The error sum of squares also increases, which is inevitable because there is one less free parameter in the model, but only does so by 0.6%. It follows that even a one-parameter model estimatingĈ Pα but setting C D0 = 0 can give a reasonably close fit to the free-flight data.

(d) Goodness of fit of the stroke-averaged aerodynamic forces
Because the blade element model is effectively fitted as a regression model forced through the origin (Eq. 2.20), its R 2 statistic is not well defined. Hence, to assess the goodness of fit of the blade element model with respect to the forces measured in each of the body axes, we regressed the fitted stroke-averaged aerodynamic forces on the measured stroke-averaged aerodynamic forces, without forcing the regression line through the origin (Fig. 7). We did this separately for each body axis, and found that a large proportion of the variation in the measured forces was explained in both the x b -and the z b -axes (R 2 = 82.5% and R 2 = 79.0%, respectively), but that a much smaller proportion of the variation in the measured forces was explained in the y b -axis (R 2 = 18.0%). The regression slope was close to one in the x b -and the z b -axes (slopes of 0.92 and 1.05, respectively), but was greatly attenuated by the noise in the y b -axis (slope of 0.62). In contrast, the regression intercept was close to zero in all three body axes (x b : −0.044 mN; y b : 0.000 mN; z b : −0.021 mN). The comparatively poor fit of the blade element model to the forces measured in the y b -axis reflects the fact that the total range of the lateral stroke-averaged aerodynamic forces was small (Fig. 3e), which in turn leads to a lower signal to noise ratio in the   (Fig. 7A,C). Visual inspection of the regression plots shows that the aerodynamic forces measured in the x b -axis are fitted closely over their entire range (Fig. 7A), but that the blade element model systematically under-predicts the magnitude of the largest aerodynamic forces produced in the z b -axis (bottom left region of Fig. 7C). Even so, the blade element model does a good job of fitting the time-history of the stroke-averaged forces estimated from the body's acceleration in all three axes (see representative sequence plotted in Fig. 8).

(e) 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 statistical reason to assume that the resulting model will perform well in predicting the time-varying aerodynamic forces through the wingbeat, although the logic of the underlying aerodynamic model is such that it should. On this basis, we used the blade element model to predict how the aerodynamic forces are expected to vary through a standard wingbeat in hovering flight (Fig. 9, 10; Supplementary Videos S2, S3), to allow us to assess the relative contributions of lift, drag, and added mass at different stages of the wingbeat. We determined the kinematics of this standard wingbeat as the mean through time of a set of wingbeats chosen to be representative of hovering flight. We first selected the subset of wingbeats representing equilibrium flight, defined as those wingbeats for which the vertical force was within ±5% of body weight, and for which the horizontal acceleration was < 0.5ms −2 . From this subset of equilibrium wingbeats, we then sub-selected the 1st percentile with respect to body speed, and averaged these to determine the standard hovering wingbeat (Supplementary Data S1). The final  subset of n = 265 wingbeats used for this averaging drew on data from a total of 51 different flight recordings from a total of 19 individuals, and their mean should therefore be broadly representative of hovering flight (Fig. 5). As a final check on the robustness of our predictions to errors in parameter estimation, we also 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 Section (c). Despite the variation in the aerodynamic force coefficient parameters fitted in the subsampling analysis (Fig. 6B), the resulting variation in the predicted time-varying aerodynamic forces was slight in comparison to their variation through the wingbeat (Fig. 9). 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.6% of the stroke-averaged vertical force (Fig. 9D. The added mass and drag forces each make a small net positive contribution to weight support, providing 13.8% and 8.6%, respectively, of the predicted stroke-averaged vertical force. The predicted lift force peaks at the middle of each halfstroke, when the wing's translational velocity is highest (Fig. 10F,H; Supplementary Video S3), so this is also the phase of the stroke during which the majority of the vertical force is expected to be produced ( Fig. 9D; Fig. 10B,D). In contrast, the added mass force peaks at the beginning of each half-stroke, when the wing's acceleration normal to its chord is highest (Fig. 10E,G; Supplementary Video S3). The drag force has a more complicated time-history again, peaking at several points through the stroke (Fig. 10E-H; Supplementary Video S3). Aerodynamic force production in the transverse y b -axis is qualitatively similar on both the upstroke and downstroke (Fig. 9B), and the same is true of the vertical component of the aerodynamic force (Fig. 9D), but in each case the amplitude of the forces is somewhat diminished on the upstroke relative to the downstroke. In contrast, aerodynamic force production along the x b -and z b -axes of the body displays a marked asymmetry between the upstroke and the downstroke (Fig. 9A,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 (Fig. 9D).

Discussion
The key contribution of this paper is to derive a physics-based blade element model with just two empirically fitted coefficients: an aerodynamic force coefficient derivative (Ĉ Pα ), and a drag coefficient offset (Ĉ D0 ). We estimated these two parameters for the largest available free-flight dataset of insect wing and body kinematics to date [31], comprising N = 26, 541 wingbeats from hovering and manoeuvring hoverflies. We found that the model captured approximately 80% of the variation in the measured stroke-averaged aerodynamic forces in the sagittal plane of the body. Importantly, the model takes full account of the torsional deformation of the wing, without limiting the motion to a fixed stroke plane, or treating the wing as a flat plate. This matters, because wing flexion is a characteristic feature of insect flight, which can reduce aerodynamic power requirements and enhance lift forces [20,33]. Another important feature of the 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. 14,25], 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 ignore the body's velocity and acceleration when fitting the blade element model, then the error sum of squares increases by 36%, and the estimate of the drag coefficient offset becomes unreasonably high (Ĉ D0 = 0.78, as compared toĈ D0 = 0.17 with the body motion included).

(a) 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 (Fig. 2), there are no modelling assumptions that are specific either to hoverflies or to the dataset we used. Indeed, one of the great strengths of our approach is that the same blade element model can be straightforwardly fitted to any similar data set, using linear least squares to optimise the force coefficient parameters: MATLAB code implementing the blade element model is provided as Supporting Data S1 to this end. Moreover, 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. As with any form of regression modelling, however, an important caveat is that any such dataset must contain a large enough range of variation in the independent variables to ensure a good fit to the data (compare, Fig. 7a,c with Fig. 7b). In principle, our regression estimates of the lift and drag coefficients could be replaced by prior estimates from either model wings [6] or computational fluid dynamics modelling [21]. 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 [19], may differ markedly from the aerodynamic properties of a real wing in vivo. These problems are avoided by our approach of fitting the parameters of the blade element model empirically to free-flight data from live insects.

(b) Comparison to other models
Blade element models have been used previously to predict the aerodynamic forces of insect flight, but have not usually been paired with detailed kinematic data from free-flying insects. Our model differs from those used previously in respect of: (i) the fact that it incorporates the full threedimensional motion of the wing, including torsional deformation; (ii) its method of accounting for the rotational force by calculating the angle of attack at the three-quarter chord point, as proposed by Ellington [11]; (iii) the fact that the resultant quasi-steady aerodynamic force is not constrained to act either perpendicular to the flow or normal to the wing's surface; and (iv) the simplicity of its form, developed from first principles. The first of these distinguishing features relates mainly to the availability of data, but the others are more fundamental.
Concerning our treatment of the rotational forces, most recent blade element models explicitly include a rotational lift component, sometimes called a rotational circulation force [e.g. 24,26]. Moreover, the blade element model of Nakata et al. [21], with coefficients fitted using computational fluid dynamics, also includes a rotational drag force. These rotational forces are intended to capture the aerodynamic effects of pitching or twisting about the spanwise axis of a flapping wing. However, by defining the velocity of each blade element at its three-quarter chord point as Ellington [11] originally proposed, 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 numerical terms that need to be determined [27].
Concerning the third distinguishing feature of our model, the instantaneous direction of the quasi-steady pressure force is determined 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 [15,23,24,26]. This approach is inconsistent with our assumption of a drag offset, however, because it implies that a circulatory force will be produced perpendicular to the chord even when the aerodynamic angle of attack is zero. Instead, the presence of a drag offset in our model implies that the direction of the quasi-steady aerodynamic force must vary with respect to the chord. We can quantify the effect of this drag offset by examining how the angle γ between the predicted quasi-steady aerodynamic force and the local chord changes with angle of attack (Fig. 11A). At angles of attack above about 45 • , the quasi-steady force approaches a 90 • angle to the chord. At lower angles of attack, the quasi-steady force vector is tilted back, becoming tangent to the blade element at zero angle of attack. 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 idealised thin aerofoil, they do in fact have a finite thickness, particularly at the leading edge, which is typically structured as a reinforced spar to reduce unwanted bending and torsion. Wing corrugation due to venation also makes wings inherently three-dimensional structures, which can further increase profile drag [28,30]. Aerodynamic measurements of real insect wings [9,10], or mechanical and computational models thereof [6,22], 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. Whilst it is true that insect wings 20 operate at characteristically high angles of attack, in the hoverfly data that we have presented here, 63% of all of the aerodynamic force produced by the wing is predicted to occur at angles of attack α < 45 deg (Fig. 11B), with a mode at α = 33 deg, for which the force vector will be tilted back significantly from 90 deg (Fig. 11A).
Perhaps the most important feature of our model, though, is the simplicity of its form, which was developed from first principles, rather than heuristically from the data. Remarkably, this approach leads to a model with only two fitted parameters, and in turn provides physical insight into the more complicated form of other models that have been fitted previously. In particular, the empirical model fitted by Dickinson et al. [6] to force measurements from a robotic flapper modelled the lift and drag coefficients as C L = 0.225 + 1.58 sin (2.13α − 0.14) and C D = 1.92 − 1.55 cos (2.04α − 0.17), where α is in radians. These formulae look odd mathematically, but they make more sense physically when it is seen that they approximate trigonometric double angle formulae. Specifically, noting that sin (2α) = 2 sin α cos α and cos (2α) = 1 − 2 sin 2 α, we can rewrite Eqs. 3.1 and 3.2 of our own fitted model as: These equations are quite similar to those fitted by Dickinson et al. [6], from which we conclude: (i) that the equations fitted by Dickinson et al. [6] approximate trigonometric double angle formulae; and (ii) that they do so because 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. We conclude that the model by Dickinson et al. [6] is likely to have been overfitted.

(c) Limitations
As noted in the Results, the blade element model systematically under-predicts the magnitude of the largest aerodynamic forces produced in the z b -axis (bottom left region of Fig. 7C). This discrepancy presumably indicates a nonlinearity in the aerodynamic forces at high angle of attack that the quasi-steady blade element model fails to capture. Most of the aerodynamic force in the z b -axis is produced during the downstroke (Fig. 9B). We therefore hypothesise that an unsteady mechanism associated with the downstroke that is not captured by the blade element model may become increasingly important as insects produce stronger flight forces. Two candidate unsteady mechanisms include wake capture from the previous wingbeat, and wing-wing interactions in the form of clap-and-fling [11]. Both effects may be enhanced to enable high force production, but whereas wing-wake interactions would be expected to be present on both the downstroke and upstroke, wing-wing interactions will only affect the start of the downstroke. This is because the wingtips only approach each other towards the end of the upstroke, and do not approach one another closely at the end of the downstroke (Fig. 12A). Interestingly, although there is a positive association between the stroke amplitude and the magnitude of the measured aerodynamic force along both the x b -and z b -axes (Fig. 12B), as would be expected under a quasi-steady model of the forces, this association is much stronger in x b (R 2 = 0.53) than in z b (R 2 = 0.27). This suggests that increases in the magnitude of the aerodynamic forces in z b 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 (Fig. 12C), which is stronger in the z baxis (R 2 = 0.62) than in the x b -axis (R 2 = 0.25). This negative association would be expected under an unsteady clap-and-fling mechanism, and the strength of the association in z b 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 force production in the z b -axis is affected more by wing-wing interactions than force production in the 21 Figure 12. Histograms of wing tip separation and stroke amplitude, and their association with the measured strokeaveraged 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 downstoke, versus the measured forces in the x b -(blue) and z b -(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, normalised by the mean wing chord; R 2 statistics for the linear associations between the wingbeat parameters and the measured forces are shown for each axis.
x b -axis. Since the blade element model will not capture this nonlinearity directly, but will instead attempt to force a straight line through the data, it follows that the linear parameter estimate for 22 C Pα can be expected to overestimate the true value of the quasi-steady force coefficient derivative C Pα .
Finally, although the effects of body motion are captured in our modelling of the wing kinematics and aerodynamics, the body itself will produce drag, and perhaps some lift, in forward flight [10]. 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, the aerodynamic forces on the body are likely to have been overwhelmed by the aerodynamic forces acting on the wings.

Conclusions
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 freeflying insects in manoeuvring flight. The scope of this work extends beyond proof-of-principle, however. The alternative approach of 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. 4]. However, it is also 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. 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 and transparent physical basis of the underlying aerodynamic equations; (ii) the complexity of the wing kinematics and body motions that it models; and (iii) the fact that the aerodynamic force coefficients are fitted empirically to freeflight data from real insects, implicitly capturing the full scope of the insect's aerodynamic force production. We therefore expect that it will serve as a useful, validated tool for future research on insect flight dynamics and control.
Data Accessibility. The supplementary data and code supporting this article will be made available through figshare.
Authors' Contributions. SW collected the published dataset on which this analysis is based, wrote code, analysed data, prepared figures, and co-wrote the paper. GT conceived the study with SW, wrote theory, analysed data, and co-wrote the paper. Both authors gave final approval for publication and agree to be held accountable for the work performed therein.
Competing Interests. We have no competing interests. angle between local velocity of wing and its surface γ angle between the predicted quasi-steady aerodynamic force and the local chord θ deviation angle, corresponding to elevation of wing tip in body-fixed axes φ stroke angle, corresponding to azimuth of wing tip in body-fixed axes ω 0 wing pitch angle offset ωr linear twist gradient of wing ω(r) local pitch angle of wing at radial coordinate r Forces A elemental added mass force (vector) D, D elemental quasi-steady drag force (vector, scalar) F total friction force (vector) L, L elemental quasi-steady lift force (vector, scalar) P total pressure force (vector) Pqs elemental quasi-steady pressure force (scalar) Pus elemental unsteady pressure force, or added mass force (scalar) Axes {x b , y b , z b } right-handed, body-fixed axis system with origin at body centre of volume {x L , y L , z L } right-handed, body-fixed axis system with origin at base of left wing {x R , y R , z R } right-handed, body-fixed axis system with origin at base of right wing {X, Y, Z} right-handed, lab-fixed axis system 1 ⊥S unit vector normal to blade element surface, resolved in the body axes 1 U unit vector parallel to velocity of blade element, resolved in the body axes 1 ⊥U,r unit vector perpendicular to blade element velocity and span, resolved in the body axes 26 Other symbols c wing chord length c 3/4 , c 1/2 distance of point back from spanwise axis (three-quarter chord point, half-chord point) C D quasi-steady drag coefficient C D0 quasi-steady drag coefficient offset C L quasi-steady lift coefficient C P quasi-steady pressure force coefficient C Pα quasi-steady pressure force coefficient derivative g gravitational acceleration in an inertial frame of reference i blade element number m insect mass n wingbeat number N total number of wingbeats Q ⊥S local flow speed normal to surface in an inertial frame of reference r radial coordinate along wing span R rotation matrix bringing vector from the body axes into the lab axes