SUMMARY
Recent studies suggest that fruit flies use subtle changes to their wing motion to actively generate forces during aerial maneuvers. In addition, it has been estimated that the passive rotational damping caused by the flapping wings of an insect is around two orders of magnitude greater than that for the body alone. At present, however, the relationships between the active regulation of wing kinematics, passive damping produced by the flapping wings and the overall trajectory of the animal are still poorly understood. In this study, we use a dynamically scaled robotic model equipped with a torque feedback mechanism to study the dynamics of yaw turns in the fruit fly Drosophila melanogaster. Four plausible mechanisms for the active generation of yaw torque are examined. The mechanisms deform the wing kinematics of hovering in order to introduce asymmetry that results in the active production of yaw torque by the flapping wings. The results demonstrate that the strokeaveraged yaw torque is well approximated by a model that is linear with respect to both the yaw velocity and the magnitude of the kinematic deformations. Dynamic measurements, in which the yaw torque produced by the flapping wings was used in realtime to determine the rotation of the robot, suggest that a firstorder linear model with strokeaverage coefficients accurately captures the yaw dynamics of the system. Finally, an analysis of the strokeaverage dynamics suggests that both damping and inertia will be important factors during rapid body saccades of a fruit fly.
INTRODUCTION
Flying insects use a variety of mechanisms (Sane, 2003) in order to control forces and moments during flight, many of which are likely to involve subtle changes in wing kinematics. While our understanding of the basic mechanisms behind aerodynamic force production has increased considerably in recent years (WeisFogh, 1973; Ellington et al., 1996; Dickinson and Götz, 1996; Dickinson et al., 1999; Sane and Dickinson, 2001), there still exists a substantial gap in our understanding of how changes in wing kinematics translate into changes in forces and moments. In addition, understanding the relationship between wing kinematics and the production of forces and moments is only a first step, as ideally we would like to know how animals actively regulate wing motion in order to control their overall trajectory through space. This problem is complicated by the fact that the forces and moments produced by flapping wings in turn depend upon the motion of the animal, i.e. upon its rotational and translational velocities. Thus, in order to tackle this problem, we require a detailed understanding of how forces and moments depend upon changes in both wing kinematics and body motion.
Dynamically scaled robotic models have been a useful tool in the investigation of flapping flight aerodynamics (Usherwood and Ellington, 2002a; Usherwood and Ellington, 2002b; Dickinson et al., 1999; Sane and Dickinson, 2001). They have allowed researchers to examine many of the aerodynamic mechanisms involved in force production such as the leading edge vortex, wakecapture and clapandfling (Maxworthy, 1979; Ellington et al., 1996; Birch and Dickinson, 2003). In addition, robotic models have enabled researchers to systematically examine the effects of kinematic variables such as the timing of wing rotation, stroke amplitude or wing–wing interactions (Sane and Dickinson, 2001; Sane and Dickinson, 2002; Lehmann et al., 2005). In all these prior examples, the robot is used as a simple feedforward device, in which kinematics are prescribed and the resulting forces and flows are measured. Another class of robotic systems, called ‘captive trajectory systems’, is used by aerodynamics researchers to investigate the freeflight dynamics of vehicles and to determine the trajectory of objects separated from a vehicle during flight (Woods et al., 2001; Ahmadi et al., 2005; Guigue et al., 2006; Guigue et al., 2007). A captive trajectory system consists of a robotic mechanism for moving the model, sensors for measuring aerodynamic forces and moments acting on the model, and a feedback mechanism for moving the model via the robotic mechanism. The motion of the model is determined in realtime, based on its mass and moments of inertia, using the measured forces and moments.
The wing and body kinematics of freeflying fruit flies performing rapid stereotyped turns called body saccades – during which flies change orientation by 90 deg in less than 50 ms – have been measured using highspeed videography (Fry et al., 2003; Fry et al., 2005). These prior studies suggest that changes in yaw torque were produced mainly by subtle changes in wing stroke amplitude and adjustments of the stroke plane. Recently, automated methods of extracting wing and body kinematics from multicamera highspeed video sequences have been developed (Fontaine et al., 2009; Ristroph et al., 2009), which should enable researchers to process much larger data sets and provide further insight into exactly how fruit flies modulate their wing kinematics to control forces and moments. Such methods have recently been used to study how flies modulate the pitch of their wings in order to induce sharp turns (Bergou et al., 2009). However, at present, the complete repertoire of kinematic changes used to control forces and moments during body saccades or other maneuvers is not known.
Early studies of yaw turns in Drosophila (Reichardt and Poggio, 1976) suggested that dynamics of the turns were dominated by friction. However, a subsequent study of freeflight saccadic turns in Drosophila (Fry et al., 2003) came to a different conclusion, suggesting that inertia played a dominant role during rapid flight maneuvers. This conclusion was based on an analysis using a yaw damping coefficient estimated by integration of Stoke's law for the fly body only. Also, the additional damping caused by the flapping wings themselves was ignored due to the fact that the aerodynamic forces were measured using wing motion in bodycentered coordinates. Recently, studies using computational fluid dynamics to examine the time course of torque production during freeflight saccades (Ramamurti and Sandberg, 2007) have raised questions in regard to this conclusion. An estimate of the yaw damping due to the flapping wings, based on a quasisteady analysis (Hesselberg and Lehmann, 2007), suggested that the damping due to the flapping wings was approximately two orders of magnitude larger than that for the body, implying that yaw turn dynamics are dominated by damping rather than inertia. A later, but similar, quasisteady estimate (Hedrick et al., 2009) also suggested that the yaw damping due to the flapping wings is significantly larger than that for the body. In both these quasisteady estimates, it was shown that the yaw damping due to the flapping wings should be proportional to the yaw velocity, rather than being proportional to the square of the yaw velocity as might be expected by a simpler model. The rotational damping due to the flapping wings was recently investigated experimentally (Cheng et al., 2009) using a dynamically scaled robotic wing. The authors found a substantial difference in the yaw torque produced by the flapping wing in the presence and absence of body rotation. In addition, these authors used a quasisteady model to estimate the damping coefficients about the roll, pitch and yaw axes as well as cross coupling between the axes. The results suggested that the damping torque is most prominent about the turning axis upon which the angular velocity is specified and that the damping coefficients for the roll and yaw axes are roughly twice as large as that for the pitch axis.
In the present study, we examine the dynamics of yaw turns during hovering conditions using a dynamically scaled robotic model that utilizes active feedback linking torque production to body motion. We consider four biologically plausible mechanisms for actively generating yaw torque that take the form of deformation modes that introduce kinematic asymmetry into a set of baseline wing kinematics. The magnitude of the kinematic asymmetry for each mode is controllable via a single deformation parameter. The yaw torque, for each mode, is measured as a function of yaw velocity and the magnitude of the deformation parameter. In addition, the robotic model is equipped with a torque feedback mechanism that is essentially a single degreeoffreedom captive trajectory system for motion about the yaw axis. The torque feedback mechanism enables the yaw torque produced by the robotic model to be measured and used in real time to set the yaw velocity of the model by simulating the inertial dynamics of the fly. Using the torque feedback system, the yaw turn dynamics of a fruit fly were estimated by examining input–output behavior of the system in response to sinusoidal and squarewave inputs.
Our measurements show that the strokeaveraged yaw torque produced by the flapping wings is about 68 times larger than that for the body, in agreement with the quasisteady estimates in previous studies (Hesselberg and Lehmann, 2007; Hedrick et al., 2009; Cheng et al., 2009). The strokeaveraged torque is found to be linear with respect to yaw velocity, and the slope of this dependence is found to be independent of the choice of mechanism for active torque generation. In addition, the active strokeaveraged yaw torque generated by the four deformation modes is found to depend linearly on both the deformation parameters and the yaw velocity. The dynamic measurements, using the torque feedback mechanism, demonstrate that a firstorder linear strokeaveraged model captures the yaw turn dynamics of a fruit fly with reasonable accuracy. Finally, the dynamic tests suggest that, despite the relatively higher damping coefficient due to the flapping wings, both inertia and damping play a significant role in the dynamics of saccadic turns.
MATERIALS AND METHODS
Robotic fly apparatus
We designed a flapping twowinged robotic apparatus, similar to that described previously (Dickinson et al., 1999), in which the entire wing assembly was capable of rotation about the yaw axis (Fig. 1). The mechanism was designed to enable quantitative measurements of the yaw torque produced by wings that were simultaneously flapping and rotating about the yaw axis. In addition, the robotic mechanism provided realtime feedback capability in which the torque produced by the wings was used to actuate the robot via a torquefeedback mechanism. When operated in this mode, the torque produced by the wings, measured using a torque sensor, was used to specify the yaw rotation rate in real time using a dynamic model. The torquefeedback mechanism provides several key advantages over allowing the system to freely rotate about the yaw axis on bearings. First, it enables us to match the dimensionless yaw moment of inertia of the fly without resorting to a large external mass. Second, it allows the yaw moment of inertia to be easily changed in software. Third, it makes it possible to emulate a truly frictionless bearing for motions about the yaw axis. Finally, it enables the use of the same mechanism for both prescribed yaw motions and for torque feedback.
The robotic mechanism consisted of a wing assembly connected to a yaw rotation stage by a 0.21 m shaft. The shaft was mounted in a pillow block bearing on the yaw rotation stage and connected to the wing assembly through a torque sensor mounted on the base plate of the assembly. The motion of the yaw rotation stage was controlled by a stepper motor (M22183.0S; Schneider Electric Motion, Marlborough, CT, USA) that was connected to the shaft via a timing belt. The wing assembly contained two wing mechanisms, each of which consisted of an array of three motors mounted on the base plate of the assembly. The degrees of freedom for each wing were controlled independently by one motor of the array and are illustrated in Fig. 1. The stroke position of each wing was controlled by a stepper motor (M17151.5D; Schneider Electric Motion) whereas the rotation and deviation of each wing were controlled by RC servo motors (HSC5996TG; Hitec RCD, Poway, CA, USA).
The wings were immersed in a 1 m×2.4 m×1.2 m tank filled with mineral oil (Cheveron Superla white oil; Chevron Texaco Corp., San Ramon, CA, USA) of density 880 kg m^{–3} and kinematic viscosity 115 cSt at room temperature. Custom software written in Python and C enabled control of the robot from a PC. A torque sensor (TQ20225Z; Omega Engineering, Stamford, CT, USA) mounted on the base plate of the wing assembly measured the yaw torque produced by the rotating and flapping wings. The torque sensor had a full scale range of 0.175 Nm and an accuracy of 0.2% fullscale output. The isometrically enlarged wings of the robotic model were based on the planform of a Drosophila melanogaster wing. The wings of the robotic model were cut from an acrylic sheet and had the following dimensions: length (R)=0.23 m, mean chord , width=0.0023 m. The rotation axes of the two wings were parallel and separated by 0.11 m.
Baseline wing kinematics
The baseline wing kinematics used in this study are based on the treatment given previously (Berman and Wang, 2007) and were chosen to be an idealized representation of the wing kinematics of D. melanogaster. The stroke position (ϕ_{b}), deviation (θ_{b}) and rotation angle (α_{b}) for the baseline kinematics are given as follows: (1) (2) (3) where f is the flapping frequency, ϕ_{0} is the stroke amplitude, α_{0} is the rotation amplitude, and the parameters k_{ϕ} and k_{α} control the shape of the kinematics. The parameter k_{ϕ} provides an experimental variable that can distort the shape of the stroke position waveform from a sinusoid to a triangle waveform as k_{ϕ} varies from 0 to 1. Similarly, the parameter k_{α} provides an experimental variable that can distort the rotation waveform from a sinusoid to a step function as k_{α} varies from 0 to ∞. Values of k_{ϕ}=0.01 and k_{α}=1.5 were selected to produce waveforms that resemble an idealized version of the wing kinematics of D. melanogaster. Similarly, a value of ϕ_{0}=70 deg was used to give a peaktopeak stroke amplitude of 140 deg, and a value of α_{0}=45 deg was used to give a 45 deg angle of attack at midstroke.
Kinematic deformation modes
The kinematic deformation modes modify the baseline wing kinematics by introducing asymmetry between the left and right wings. The deformation modes were selected as plausible mechanisms for generating yaw torque during flapping flight. Four deformation modes were considered: differential angle of attack, differential deviation, differential stroke plane rotation and differential stroke velocity.
The differential angle of attack deformation mode modifies the baseline wing kinematics by introducing a difference in the angle of attack between the right and left wings that reverses on the upstroke and downstroke. This deformation mode only affects the rotation angle of the baseline wing kinematics; the stroke position and deviation angles are unaffected. The modified rotation angles for the left and right wings are given as follows: (4) where p_{a} is the deformation parameter. When p_{a} is equal to zero, the wing kinematics are symmetric; a nonzero value for p_{a} introduces asymmetry into the wing kinematics. Example left and right wing tip trajectories for kinematics deformed via the differential angle of attack deformation mode are shown in Fig. 2A.
The differential deviation deformation mode modifies the baseline wing kinematics by adding a sinusoidal variation to deviation angles of the left and right wings. In addition, the phase of sinusoid added to the deviation of the left wing is 180 deg out of phase with the sinusoid added to the deviation of the right wing. This deformation mode only affects the deviation angle; the stroke position and rotation angles are unaffected. The modified deviation angles for left and right wings are given as follows: (5) where p_{d} is the deformation parameter. Example left and right wing tip trajectories for kinematics deformed by the differential deviation deformation mode are shown in Fig. 2B.
The differential rotation deformation mode modifies the wing kinematics by differentially rotating the baseline kinematics of each wing about the lateral axis running through the two wing joints. This rotation effectively changes the stroke planes of the wings by tilting them in opposite directions. The deformation affects all three kinematic angles. Firstorder approximations of the modified angles for the left and right wings are given as follows: (6) (7) (8) where p_{r} is the deformation parameter. Example left and right wing tip trajectories for kinematics deformed by the differential rotation deformation mode are shown in Fig. 2C. The exact expressions for the modified kinematic angles are fairly complicated (see Appendix A). The approximations given in Eqns 6, 7 and 8 are accurate to within a degree for p_{r} within ±π/12 (±15 deg).
The differential velocity deformation mode modifies the wing kinematics by differentially varying the velocity of the left and right wing kinematics. This deformation mode only varies the stroke position angle. The modified stroke position angle is given as follows: (9) where (10) and p_{v} is the deformation parameter. Example left and right wing tip trajectories for kinematics deformed by the differential velocity deformation mode are shown in Fig. 2D. When p_{v} is equal to zero, the wing kinematics are symmetric. When p_{v} is greater(less) than zero, the left wing moves faster(slower) on the downstroke and slower(faster) on the upstroke whereas the right wing moves slower(faster) on the downstroke and faster(slower) on the upstroke.
In experiments in which the four kinematics deformation modes were combined, the differential angle of attack, deviation and velocity modes were applied first, followed by the differential rotation mode. Complete expressions for the wing kinematics deformed using the combinations of the deformation modes are given in Appendix B.
Measurements of strokeaveraged yaw torque
Measurements of the strokeaveraged yaw torque for each kinematic deformation mode were made as a function of the yaw velocity and deformation parameter. During these measurements, the robotic mechanism was set to rotate at a constant yaw velocity, and the value of the deformation parameter was held constant throughout the trial. Each trial consisted of five flapping cycles, and the yaw torque was averaged over the middle three cycles, as shown in Fig. 3. The trials were conducted at yaw velocities of –7 to 7 deg s^{–1} in steps of 1 deg s^{–1} for the values of deformation parameters summarized in Table 1.
Two dimensionless parameters are required in order to achieve an accurate dynamic scaling of the torques obtained via the robotic model: the Reynolds number (Re) and the dimensionless yaw velocity (ω*). The Reynolds number is given by: (11) and the dimensionless yaw velocity is given by: (12) where R is the wing length, Φ is the (peaktopeak) stroke amplitude, f is the flapping frequency, is the mean wing chord, and ω is the yaw velocity. A flapping frequency of 0.167 Hz was selected for the wing kinematics to give a Re of approximately 100, matching the value appropriate for D. melanogaster (Lehmann and Dickinson, 1997). The range of yaw velocities of the trials corresponds to yaw velocities of –8400 to 8400 deg s^{–1} for a fruit fly and spans the average yaw velocity for a saccade of approximately 2000 deg s^{–1} (Tammero and Dickinson, 2002; Fry et al., 2003).
The measured yaw torque can be written in dimensionless form as: (13) where ρ is density. Using the fact that the dimensionless yaw torque for the robotic model and for the fruit fly must be the same, the yaw torque measured by the robotic model is related to that acting on a fly according to: (14) where the subscript ‘fly’ and ‘robo’ refer to variables corresponding to a fruit fly or the robotic apparatus, respectively. The scaling factor, in Eqn 14, relating the torque measured by the robotic model to those acting on a fruit fly is estimated to be 4.54×10^{–7}. Note, the torque relationship given in Eqn 14 does not rely on any assumptions regarding the magnitude of Re.
Dynamic measurements using torquefeedback
Measurements using the robotic model's torquefeedback mechanism were made to measure the dynamic response of a fly to the torque produced via asymmetry in the wing kinematics. When in torquefeedback mode, the robotic model uses the instantaneous measured yaw torque to determine the heading and yaw velocity in real time by integrating a dynamic model of the fly's inertial dynamics. The model of the fly's inertial dynamics is given by: (15) (16) where ψ is the heading angle, I is the moment of inertia about the yaw axis, ω is the yaw velocity, τ_{meas} is the yaw torque measured by the sensor, and b is an optional damping term. This equation was integrated using the classical RungeKutta method (Butcher, 2003) to set the yaw velocity and heading angle of the system at each time step of the realtime loop, which was updated at 3 kHz. For all of the experiments described in this manuscript, the optional damping term, b, was set to zero, so that all of the damping came from the aerodynamic forces acting on the wings and body of the robotic model.
Two dimensionless parameters are required in order to obtain an accurate dynamic scaling of the torque and yaw velocity for the dynamic measurements: Re, as defined in Eqn 11, and the dimensionless moment of inertia about the yaw axis, I*. The dimensionless moment of inertia is given by: (17) where ρ is the density of the fluid.
As with the strokeaveraged torque measurements, all experiments were performed with a flapping frequency of 0.167 Hz to give a Re of 100. Similar to Fry et al. (Fry et al., 2003), the moment of inertia of a fly about the yaw axis was modeled as a cylinder with length 2.44×10^{–3} m, radius 4.28×10^{–4} m, and a mass of 1.1×10^{–6} kg (Lehmann and Dickinson, 1997). The body posture of a Drosophila in hovering flight was mimicked by orienting the axis of the cylinder at an angle of 35 degrees with respect to vertical. The moment of inertia about a vertical axis running through the center of mass of the cylinder was calculated to be 2.77×10^{–13} kg m^{–2}, and the dimensionless moment of inertia is given by 1.97×10^{3}.
The yaw torque measured by the robotic model can be related to the yaw torque acting on a fruit fly using Eqn 14. The relationship between the yaw velocity of the robotic model and the yaw velocity of a fruit fly is given by: (18) Similarly, the relationship between time for the robotic model and time for the fruit fly is given by: (19) The scaling factor in Eqns 18 and 19 (f_{robo}/f_{fly}) is equal to 1200.
Two types of measurements were made to determine the dynamic response of the fly to yaw torque produced by asymmetry in the wing kinematics. In the first type of measurement, for each deformation mode, the deformation parameter was varied in a sinusoidal fashion as a function of time at a fixed amplitude and frequency. The amplitudes of the sinusoidal variations were scaled so that each deformation mode produced approximately the same strokeaveraged yaw torque as a function of time and yaw velocity as follows: (20) where x indicates the deformation mode, s_{x} is the scaling factor for deformation mode x, and u_{cos}(t) is the sinusoidal variation. The function u(t) is given by: (21) where u_{0} is the amplitude of the waveform, and f_{u} is the frequency of the waveform. The appropriate values for the scaling factors were determined from the results of the strokeaveraged torque measurements using a linear regression and are summarized in Table 2.
The sinusoidal trials were performed with an amplitude u_{0}=0.12 at 25 frequencies, ranging from f_{u}=f/2 to f_{u}=f/60 with logarithmic spacing. Thus, in dimensionless form: (22) where the frequencies ranged from 1/2 to 1/60. Trials consisted of five cycles of the sinusoidal waveform. Example data for a single cycle of a representative trial are shown in Fig. 4.
In the second type of measurement, for each deformation mode, the deformation parameter was a squarewave with a fixed amplitude and frequency. As with the sinusoidal measurements, the amplitudes of the squarewaves (u_{sqr}) were scaled so that each deformation mode produced approximately the same strokeaveraged torque as a function of time and yaw velocity: (23) where s_{x} represents the scaling factors given in Table 2, and u_{sqr}(t) is the squarewave function. In order to avoid discontinuities in the wing kinematics, an approximation of a squarewave, which transitions linearly from hightolow and lowtohigh over a short time period, was used. A complete expression of the approximate squarewave function is given in Appendix C. Squarewave amplitudes of 0.05 and 0.1 were selected to produce turns with approximately 1 and 2 times the average velocity of a saccade, i.e. the first amplitude produces a roughly 90 deg turn in 10 wing strokes whereas the second produces a roughly 180 deg turn in 10 wing strokes. By the same reasoning, the frequency of the squarewave was set to 120 s or 20 wing strokes. Thus, the robotic model produced positive torque for 10 wing strokes and negative torque for another 10 wing strokes. Trials consisted of five cycles of the squarewave. Example data for a single cycle of a representative trial are shown in Fig. 5.
For both types of dynamic measurements, sinusoidal and squarewave, the time course of the resulting yaw velocities was used to access the dynamic response of the fly to torque produced by asymmetries in the wing kinematics. For the sinusoidal variations in the deformation parameters, the resulting yaw velocities were approximately sine waves. For each input frequency, the amplitude and phase (relative to the input sine wave) of the resulting yaw velocities were estimated using a leastsquares fit. Using the amplitudes, the gain (in dB) of the system as a function of frequency was calculated as: (24) where is the amplitude of the dimensionless yaw velocity ω*.
For the squarewave measurements, the resulting yaw velocities are compared with the predictions of a firstorder linear dynamic model of the system dynamics based on strokeaverage coefficients: (25) where and are the strokeaveraged dimensionless damping and actuation coefficients, respectively. The dimensionless damping coefficient is given by: (26) and the dimensionless actuation coefficient is given by: (27)
Measurement of body damping
In order to determine the relative magnitude of the rotational damping (countertorque) due to the fly's body compared to that produced by the flapping wings, measurements were made of the rotational damping on a model fly body (Fig. 6) and on the robotic model with the wings removed. The body model was constructed using SolidWorks (Dassault Systèmes SolidWorks Corp., Concord, MA, USA), from images of female Drosophila melanogaster, and was fabricated in ABS plastic using fused deposition modelling. The length of the model, measured along the long body axis, was 0.27 m.
The countertorque produced by the body model and the robotic model (without wings) was much smaller than that due to the flapping wings, and a higher viscosity oil was required in order to obtain large enough torques for accurate measurements. The experiments were performed in mineral oil (Cheveron Superla white oil; Chevron Texaco Corp.) with a kinematic viscosity of 1.87 St and a density of 860 kg m^{–3}. For the experiments using the body model, the wing assembly was removed from the rotation stage and replaced with the body model. In each trial, the body or robotic model (without wings) was set to rotate at a constant yaw velocity from a heading of 0 to 360 deg. The yaw torque was then averaged from the segment of data from 180 deg to 325 deg. Trials were conducted at yaw velocities of 0 to 40.7 deg s^{–1} in steps of 4.07 deg s^{–1}.
A single dimensionless parameter, Re_{b}, is required in order to obtain an accurate dynamic scaling of the torques obtained via the model fly body and the robotic model (without wings) in the higher viscosity oil. Re_{b} is the Reynolds number of the body for rotations about the yaw axis and is given by: (28) where R_{b} is the length of the body, and v is the kinematic velocity. The relationship between the yaw velocities for experiments performed at two different length scales and kinematic viscosities is given by: (29) where the subscripts 1 and 2 refer to experimental setup 1 and 2, respectively. The scaling factor in Eqn 29 relating the yaw velocities performed in the high viscosity oil to those in the low viscosity oil, is 0.073. Thus, yaw velocities for the body damping experiments, performed in the high viscosity oil, correspond with yaw velocities of 0–3 deg s^{–1} in the low viscosity oil with a step size of 0.3 deg s^{–1}. The scaling factor in Eqn 29 relating the yaw velocities performed in the high viscosity oil to those of a fruit fly is 100. Thus, the yaw velocities for the body damping experiments, performed in the high viscosity oil, correspond with yaw velocities of 0–4070 deg s^{–1} with a stepsize of 407 deg s^{–1} for a fruit fly.
The relationship between torques for experiments performed at two different length scales and kinematic viscosities is given by: (30) The scaling factor in Eqn 30 relating the torques measured in the high viscosity oil to those in the low viscosity oil is 0.0035. Similarly, the scaling factor relating torques measured in the high viscosity oil with those acting on a fruit fly is 8.1×10^{–10}. Note, the torque relationship given in Eqn 30 does not rely on any assumptions regarding the magnitude of Re.
RESULTS
Strokeaveraged yaw torque
The strokeaveraged yaw torques as a function of the yaw rate are shown, in dimensionless form, for the four deformation modes in Fig. 7. Note, for each deformation mode, the yaw torque is shown for five values of the deformation parameter for that mode. The strokeaveraged yaw torque exhibits a clear linear dependence upon the yaw rate. For each deformation mode, the offset of the linear relationship clearly depends upon the value of the deformation parameter whereas the slope of the relationship appears the same regardless of deformation mode or of the specific value of the deformation parameter. The slope of the relationship gives the strokeaveraged damping coefficient (). The mean value of the strokeaveraged damping coefficient over all trials, estimated using a leastsquares fit, was –6.4×10^{2}, and the maximum variation in the slope over all trials was less than 5%. Thus, it is clear that strokeaveraged yaw torque, for a given deformation mode x, may be approximated as a linear function of the yaw velocity, ω*, plus an unknown function, F_{x}, of the deformation parameter, p_{x}, as follows: (31)
Fig. 8 shows the strokeaveraged yaw torque, for each deformation mode, as a function of the deformation parameter for that mode. Note, for each mode, the strokeaveraged yaw torque is shown for five values of yaw velocity. The strokeaveraged yaw torque exhibits a clear linear dependence with respect to each of the four deformation parameters, and the slope of the linear relationship does not appear to depend upon the yaw velocity for any of the deformation modes considered. The slope of this relationship gives the dimensionless strokeaveraged actuation coefficient () for each deformation mode x, such that the dimensionless torque produced by that mode at yaw velocity equal to zero is given by . The mean values of the slopes for each deformation mode are summarized in Table 3. Note, the maximum variation of the slope for each mode over all trials was less than 5%. From this, it is clear that the approximate expressions for the strokeaveraged yaw torque given in Eqn 31 may be simplified to: (32)
The expressions in Eqn 32 give us an approximation for the strokeaveraged torque produced by each deformation mode as a function of deformation parameter and yaw velocity. What is not clear, however, is whether or not these expressions can be combined to give an approximation for the torque produced by kinematics that are deformed by two or more of the deformation modes simultaneously. Using Eqn 32, the expression for torque produced by the combination of two deformation modes, x and y, is given by: (33)
Fig. 9 shows the torque as a function of yaw velocity for all pairs of deformation modes for selected values of the deformation parameters, as well as the torque predicted by Eqn 33 for each of the combined pairs. The predicted torques are seen to agree extremely well with the actual torques produced by kinematics combining the different deformation modes. Thus, the strokeaveraged torque may be approximated as a linear function of the yaw velocity and deformation parameters as: (34) The torque produced by kinematics that were deformed by three or four deformation modes simultaneously is not shown, but was also well predicted by Eqn 34.
Body damping
The expression for the strokeaveraged yaw torque given in Eqn 34 includes the yaw torque produced by the body of the robotic mechanism, which may be different from that produced by a fly body. In addition, it is of interest to be able to separate and compare the torques produced by the flapping wings from those produced by the body of a fly. For these reasons, measurements were made of the yaw torque produced by the body of the robotic mechanism and of a scale model fly body. Fig. 10 shows the dimensionless yaw torque produced by both the body of the robotic mechanism (without wings) and the model fly body as a function of yaw velocity. The yaw torque produced by both the model fly body and the body of the robotic mechanism shows a clear linear dependence upon yaw velocity. The slope of this relationship gives the damping coefficient for the body. The dimensionless damping coefficients, estimated via leastsquares fit, for the model fly body and the body of the robotic mechanism are –9.3 and –4.6, respectively. Thus, we see that the yaw torque produced by the model fly body is approximately twice that produced by the body of the robotic model. However, the yaw torque produced by both the model fly body and the body of the robotic mechanism is nearly two orders of magnitude smaller than that produced by the flapping wings. The inclusion or exclusion of the yaw torque produced by either body model represents a change of less than two percent in the value of the strokeaveraged damping coefficient, . Thus, it is clear that the torque produced by the body damping is quite small compared to that produced by the flapping wings and that the effect of the body can be safely ignored in Eqn 34.
Sinewave response
The results of the previous section demonstrate that the steadystate, strokeaveraged yaw torque is well approximated by a function that is linear with respect to both yaw velocity and the deformation parameters. This suggests that the turning dynamics of a fruit fly about the yaw axis could be modeled in the strokeaveraged sense as: (35) We tested this hypothesis by performing a sinusoidal input/output analysis using the torquefeedback mechanism of the robotic model. During these tests, for each deformation mode, the input (deformation parameter) was a cosine and the output was the resulting yaw velocity of the system. Note, as described in the Materials and methods, the deformation parameters were scaled, using Eqn 20, so that each mode would produce the same strokeaveraged torque, according to Eqn 34, as a function of time. Trials were performed over a range of input frequencies in order to measure the dynamic response of the system.
Fig. 11 shows the resulting yaw velocity for all four deformation modes for a single cycle of the deformation parameter. The resulting yaw velocities show a clear sinusoidal response at the same frequency as the sinusoidal variation of the deformation parameters, which is characteristic of a linear system. The amplitude and phase of this sinusoidal component of the resulting yaw velocities depends upon the frequency of the deformation parameter sinusoid. As would be expected for a damped linear system, the amplitude decreases with increasing input frequency. Similarly, the phase shift is small for lowfrequency inputs and increases to around 90 deg as the frequency of the input is increased. The changes in amplitude as a function of frequency are summarized in Fig. 12, which shows the gain of the system in dB, calculated using Eqn 24, as a function of frequency for each deformation mode. Also shown are the predicted gains for the linear system given by Eqn 35, with coefficients determined from the strokeaveraged measurement.
For all four deformation modes, the measured gains agree quite closely with the gains predicted by the strokeaveraged linear model when the input frequency is below approximately onefifth of the flapping frequency. At higher input frequencies, the differential deviation and differential strokeplane rotation modes show slightly higher and lower gains, respectively, than predicted by the linear model. The changes in phase as a function of frequency are summarized in Fig. 13, which shows the relative phase of the output (yaw velocity) to the input (deformation parameter) signals. Also shown are the predicted phases of the linear system, given by Eqn 35, with coefficients determined from the strokeaveraged measurements. Again, the agreement of the measurements with the linear model is quite good. However, there are some differences in the behavior of the systems as compared with the linear model at higher frequencies. When using the differential angle of attack and differential rotation deformation modes, the system shows slightly less phase shift than predicted by the strokeaveraged linear model at higher frequencies. By contrast, when using the differential rotation and differential velocity deformation modes, the system shows slightly greater phase shift than predicted by the strokeaveraged linear model at higher frequencies. Overall, the results suggest that the strokeaveraged linear model provides a very reasonable approximation for the yaw turning dynamics of a fruit fly.
As evident in Fig. 11, the yaw velocity, for all four modes, shows evidence of a signal component at around twice the wing beat frequency, which is due to the unsteady modulation of the torque produced by the flapping wings. This effect is more noticeable at higher frequencies and particularly noticeable for the differential velocity deformation mode. This is due to the fact that, for this mode, the wings have unequal velocities during portions of the stroke, resulting in unsteady torques whose peak magnitude is larger than that of the other modes. Note, however, that the strokeaveraged values of the torques for this mode are approximately equal to those of the other three modes. It is interesting to investigate the origin of these higher frequency components of the resulting yaw velocities. The strokeaveraged model assumes that, at a given yaw velocity, the flapping wings will produce a constant yaw torque for a given value of the deformation parameter. This is not strictly true, as the flapping wings do not produce a steady yaw torque but rather a torque that increases and decreases throughout the stroke cycle, i.e. the torque is modulated by a function that is periodic with respect to the flapping frequency. This may be written as: (36) where g(t*) is a function that is periodic with respect to the dimensionless flapping frequency f*. The periodic function g(t*) may be written as a Fourier series: (37) where a_{n} and b_{n} are unknown coefficients. If is a strokeaveraged actuation coefficient then we must have a_{0}=2 to recover the strokeaveraged yaw torques. When the input function is a sinusoid, u(t*)=u_{cos}(t*), then: (38) The products of sinusoids in Eqn 38 may be written as the sum of sinusoids as follows: (39) and (40)
Thus, we see that the frequencies of the sinusoids in our modulated input signal are and , where n=1,2,...,∞. If the underlying yaw turning dynamics are linear, we would expect these frequency components to be present in our output signal, the yaw velocities. Fig. 14 shows the power spectral density (PSD) of the yaw velocities for a representative trial with the differential velocity deformation mode. Clear peaks are seen in the PSD at the frequencies , and . The higher frequencies, , are not shown, but the PSD also shows clear peaks at these values. The PSDs for the other deformation modes show similar peaks at these frequencies. Thus, it is clear that the additional frequencies seen in the yaw velocities are due to the modulation of the yaw torque by the flapping wings. Also, as the actuation frequency () is increased, the frequency is decreased and becomes more prominent in the output, a result that is expected for firstorder damped linear system as the gain of the system increases with decreasing frequency. Therefore, we expect the effect of the torque modulation on the output to be more pronounced for higher input frequencies, .
Squarewave response
As a second test of the yaw turning dynamics we performed an input/output analysis using square waves. During these tests, the deformation parameter or input for each deformation mode was an approximate squarewave and the output was again the yaw velocity of the system. The trials were performed with squarewave amplitudes that produced turns of approximately 1 and 2 times the average velocity of a saccade. Note that both the amplitude and frequency of the squarewave affect the average turn velocity. As with the sinusoidal input/output analysis, the deformation parameters were scaled, according to Eqn 20, so that each mode produced the same strokeaveraged torque. Fig. 15 shows the yaw velocities for all four deformation modes for a single cycle of the squarewave input in the trials that produced turns at approximately twice the average velocity of a saccade.
We also show the predictions of the squarewave response made using the strokeaveraged model given by Eqn 35. The time course of yaw velocities is characteristic of a firstorder damped linear system, and the linear strokeaveraged model with coefficients derived from the strokeaveraged measurements predicts the time course of the yaw velocities quite closely. Higher frequency oscillations, at frequencies near the flapping frequency, are apparent on the yaw velocities for all four deformation modes, reflecting the underlying modulation of the yaw torque by the flapping wings. Again, these oscillations are more prominent for the differential velocity deformation, reflecting the greater amplitude of torque modulation for this mode.
DISCUSSION
There are several possible strategies for studying how insects might modulate their wing kinematics in order to actively control forces and moments and ultimately their flight dynamics. One approach (Willmott and Ellington, 1997; Fry et al., 2003; Card and Dickinson, 2008; Liu and Sun, 2008) is to carefully record and measure the wing and body kinematics of freeflying animals. The forces and moments produced can then be analysed using simulations or dynamically scaled models. This approach yields valuable direct insight into what the animals actually do in flight to modulate the forces and moments produced. A second approach (Sane and Dickinson, 2001; Usherwood and Ellington, 2002a; Usherwood and Ellington, 2002b; Dickson and Dickinson, 2004), and the one taken in this paper, is to consider more broadly the possible changes in wing kinematics that might be used to produce a particular combination of forces or moments. This latter approach makes it possible to study the relationship between kinematics, forces and moments in a systematic manner. In the present study, we have explored four possible mechanisms for active yaw torque generation. These mechanisms, called deformation modes, function by introducing kinematic asymmetry between the left and right wings. One relatively surprising finding of this study is that, with respect to strokeaveraged yaw torque production, all of the mechanisms examined could be modeled as simple linear functions of the deformation parameters. We found that it is possible to combine the deformations, i.e. apply more than one deformation to the baseline wing kinematics, and very accurately predict the yaw torque of the combined deformation given the measurements of the yaw torque for each deformation mode in isolation. In addition, the strokeaveraged dynamics of the system, up to a scaling of the actuation coefficient, did not depend upon our choice of mechanism for yaw torque generation. Scaling the actuation coefficient will scale the gain but will not affect the phase lag of the system. These two results imply that, at least with respect to yaw torque production, all of the actuation modes are essentially equivalent and, when scaled appropriately, result in equivalent strokeaveraged yaw dynamics. These findings offer insight into the control strategies used by hovering birds and insects as well as potential design principles for the control of flapping micro air vehicles.
Strokeaveraged models
The results of this study suggest that, at least when restricted to motion about the yaw axis, a strokeaveraged model is sufficient to accurately capture the yaw turn dynamics of a fruit fly. In addition, the strokeaveraged model provides a framework for determining how much the instantaneous yaw velocity might be expected to deviate from the strokeaveraged prediction due to the fact that actively generated yaw torque is not constant but rather modulated by the flapping frequency. Differences between the prediction of the linear strokeaveraged model and the system response start to increase for control inputs with frequency greater than approximately onetenth of the flapping frequency. These differences may be important during extreme maneuvers, such as escape takeoff, where the wing kinematics change on a strokebystroke basis (Fontaine et al., 2009). It is not known whether such an approach can be extended to less restrictive motion that includes roll, pitch or translation. However, if a strokeaveraged model can be applied with sufficient accuracy to a reasonable portion of the flight envelope then this will greatly simplify the mathematical description of the motion and make questions concerning flight dynamics and control more tractable. Recent results (Bergou et al., 2009) suggest that linear strokeaveraged models can successfully capture the trends seen during turning maneuvers of freeflying insects. For flapping flight in general, the conditions under which it is reasonable to assume that a strokeaveraged model will provide an acceptable model of the flight dynamics are still unclear. Further study will be required to determine when this is the case; however, the results of the present study suggest that, at least under some conditions, the strokeaveraged approach may prove successful.
Inertia and damping
Recent studies (Hesselberg and Lehmann, 2007; Hedrick et al., 2009; Cheng et al., 2009) have suggested that the counter torque produced by the flapping wings of a fruit fly in response to rotation about the yaw axis is approximately two orders of magnitude greater than the torque produced by the body of the fly. The strokeaveraged measurements of the yaw torque produced by wings and body in the present study provide further evidence that this is indeed the case; the damping coefficient due to flapping was found to be about 68 times greater than that for the body of the fly. It has been further suggested (Hesselberg and Lehmann, 2007) that damping (or friction) is the dominant factor in the dynamics of rapid yaw turns in fruit flies. As we have shown that the turning dynamics, when restricted to motions about the yaw axis, are well modeled by a firstorder linear system, we can reexamine the roles of damping and inertia with respect to the yaw dynamics of fruit flies.
In order to examine the relative importance of damping and inertia for a firstorder linear system it is necessary to specify the input frequencies of interest. That this is the case can be seen by looking at the expressions for the gain and phase lag of the system as a function of the frequency of the input sinusoid. The gain of the system – the output amplitude over the input amplitude – is given by: (41) and the phase lag of the system is given by: (42) A system will be dominated by damping when the term is negligibly small relative to . When this is the case, the gain of the system is approximately constant and given by . In addition, the phase lag of the system is approximately zero. A system will be dominated by inertia when the damping coefficient is negligibly small relative to . In this case, the gain of the system is given by and has a strong dependence upon the frequency of the input sinusoid. The phase lag of the system, in this case, will be 90 deg. In both cases, it is the magnitude of the product of the input frequency and the moment of inertia relative to the magnitude of the damping coefficient that matters. Thus, we see that any damped firstorder linear system can be considered to be dominated by damping or inertia for a suitable choice of input frequency. In order to assess the importance of damping and inertia in the yaw turning dynamics of a fruit fly, it is necessary to select a representative frequency or range of frequencies at which to examine the system's dynamics.
The freeflight trajectories of fruit flies are characterized by straight flight segments interspersed with rapid changes in heading called body saccades (Tammero and Dickinson, 2002). During a saccade, a fruit fly rotates by approximately 90 deg in approximately 10 wing strokes (Fry et al., 2003). As an idealized representation of a saccade, we consider a sinusoid with an amplitude of 45 deg and a period of approximately 20 wing strokes. Such a sinusoid represents a continuous sequence of saccades, in which the fly turns 90 deg in one direction for 10 wing strokes and then 90 deg in the other direction for 10 wing strokes and so on. An estimate of the phase lag for a representative fruit fly in response to an input sinusoid with dimensionless frequency of 1/20 can be made using values of 1.97×10^{3} for the dimensionless moment of inertia, I*, and –6.4×10^{2} for the dimensionless damping coefficient, . The phase lag, calculated using Eqn 42, is found to approximately 44 deg. Thus, the phase lag of the system is almost exactly midway between a system dominated by damping and one dominated by inertia. In addition, the gain of the system is reduced by 3 dB from that for a system that is dominated by damping. Note, a 3 dB reduction in the gain implies that the ratio of output amplitude to input amplitude is reduced to 71% of the value it would be for a system that is dominated by damping. From this, we can conclude that neither damping nor inertia dominate the dynamic response of the system to a sinusoid representative of a saccade but rather that both inertia and damping play an important role in the response.
Fruit flies exhibit considerable variation in body mass, wing size, wing shape and other morphological parameters. In addition, the mass of an individual fruit fly might vary to a large extent depending on age, hunger status, gravidity and other factors. For this reason, it is interesting to consider how sensitive our assessment of the relative importance of inertia and damping is to changes in the moment of inertia and damping coefficient, which we can do by examining how the phase lag varies in response to changes in these parameters. The phase lag depends upon the ratio of the dimensionless moment of inertia and the dimensionless damping coefficient. This ratio is equivalent to the dimensionless time constant of the linear system. Lowering the time constant will decrease the phase lag for a given input frequency, and increasing the time constant will increase the phase lag. When the dimensionless moment of inertia and damping coefficient are as given in the previous paragraph, then the time constant is equal to 3. This can be interpreted to mean that in response to a step change in the control input u, the fly will require 3.0 wing strokes in order to reach 63% of its new terminal velocity. Varying the time constant by a factor of two in either direction gives phase lags equal to 26 deg and 63 deg, respectively, in response to a sinusoid of dimensionless frequency equal to 1/20. Both of these values for phase lag clearly fall in the intermediate regime where neither damping nor inertia dominate. We might ask at what value of time constant will the phase lag for this frequency be within 5% of the value expected of a damping or inertia dominated system, i.e. within 5% of 0 deg or 90 deg, respectively. For the phase lag to be within 5% of 0 deg, the time constant must be less than 0.25, which is approximately 8% of the value of our estimate for the time constant for a representative fruit fly. For the phase lag to be within 5% of 90 deg, the time constant must be greater than 40, which is about 1350% greater than our representative time constant. Again, it seems clear that both damping and inertia will be important factors with regard to yaw turn dynamics over a fairly wide variation in both the moment of inertia and damping coefficient.
Roll and pitch
In the present study, we measured the yaw torque produced by each deformation mode as a function of deformation parameter and yaw velocity. We did not measure the roll or pitch moments produced by the different modes and how they changed with yaw velocity. While these measurements were beyond the scope of this study, it would be very interesting to compare the roll and pitch moments for the different deformation modes. In particular, one might try and identify deformation modes that produce essentially equivalent yaw torque but different roll and/or pitch moments. Given that, with regard to yaw torque, the deformation modes could be linearly combined, modes that produce different roll/pitch moments might be combined to give independent control of yaw and roll/pitch. We are currently conducting experiments to examine this important issue of coupling.
ACKNOWLEDGEMENTS
We thank Sawyer B. Fuller and Martin Y. Peek for technical advice and assistance.
APPENDIX A. EXACT EXPRESSION FOR THE DIFFERENTIAL ROTATION DEFORMATION MODE
The exact expressions for wing kinematics modified using the differential rotation deformation mode, which differentially rotates the baseline kinematics of the two wings about the lateral axis running through the wing joints by angle p_{r}, are given as: (A1) (A2) (A3) where (A4) (A5) (A6) (A7)
APPENDIX B. WING ANGLES FOR COMBINED DEFORMATION MODES
In experiments in which the four kinematics deformation modes were combined, the differential angle of attack, deviation and velocity modes were applied first followed by the differential rotation mode. Note, that the differential angle of attack, deviation and velocity modes only affect a single wing angle each, and these angles are all distinct. Thus, these three modes are independent and their application order is irrelevant. The stroke position, deviation and rotation angles for the combined kinematics are given as: (B1) (B2) (B3)
APPENDIX C. APPROXIMATE SQUAREWAVE
In order to avoid discontinuities in the wing kinematics, an approximate squarewave function, u_{sqr}, was used in the dynamic measurements. The approximate squarewave avoids discontinuities by transitioning linearly from lowtohigh and from hightolow. The approximate squarewave function is given by: (C1) where u_{0} is the amplitude, f_{u} is the frequency and δ is the fraction of the period (1/f_{u}) spent in a given lowtohigh or hightolow transition.
FOOTNOTES

Research was supported by the Army Research Office (ARO) DAAD 19003D0004, the National Science Foundation (NSF) FIBR 0623527, and The U.S. Army Research Laboratory Micro Autonomous Systems and Technology (MAST) Collaborative Technology Alliance.
LIST OF SYMBOLS AND ABBREVIATIONS
 b
 optional damping term
 mean chord
 dimensionless strokeaveraged actuation coefficient for the differential angle of attack deformation mode
 dimensionless strokeaveraged actuation coefficient for the differential deviation deformation mode
 dimensionless strokeaveraged actuation coefficient for the differential rotation deformation mode
 dimensionless strokeaveraged actuation coefficient for the differential velocity deformation mode
 dimensionless strokeaveraged actuation coefficient for deformation mode x
 strokeaveraged dimensionless damping coefficient
 f
 flapping frequency
 f*
 dimensionless flapping frequency
 f_{u}
 waveform frequency
 dimensionless waveform frequency
 F_{x}
 unknown function of the deformation mode parameter p_{x}
 g
 periodic function
 G
 gain
 I
 moment of inertia about the yaw axis
 I*
 dimensionless moment of inertia about the yaw axis
 k_{α}
 rotation angle shape parameter
 k_{ϕ}
 stroke position shape parameter
 p_{a}
 deformation parameter for the differential angle of attack deformation mode
 p_{d}
 deformation parameter for the differential deviation deformation mode
 p_{r}
 deformation parameter for the differential rotation deformation mode
 p_{v}
 deformation parameter for the differential velocity deformation mode
 PSD
 power spectral density
 R
 wing length
 R_{b}
 body length
 Re
 Reynolds number
 Re_{b}
 Reynolds number of the body for rotations about the yaw axis
 s_{x}
 scaling factor for deformation mode x
 t
 time
 T
 wing beat period 1/f
 u_{0}
 waveform amplitude
 u_{sqr}
 squarewave amplitude
 v
 kinematic velocity
 α
 rotation angle
 α_{0}
 rotation amplitude
 α_{b}
 rotation angle of baseline kinematics
 α_{L,R}
 modified rotation angle for the left and right wings
 θ_{b}
 stroke deviation of baseline kinematics
 θ_{L,R}
 modified stroke deviation for the left and right wings
 ρ
 fluid density
 τ
 yaw torque
 τ*
 dimensionless yaw torque
 τ_{meas}
 yaw torque measured by the sensor
 ϕ_{0}
 stroke amplitude
 ϕ_{b}
 stroke position of baseline kinematics
 ϕ_{L,R}
 modified stroke position angle for the left and right wings
 Φ
 (peaktopeak) stroke amplitude
 ψ
 heading angle
 ω
 yaw velocity
 ω*
 dimensionless yaw velocity
 amplitude of dimensionless yaw velocity
 © 2010.