Dimensional analysis of spring-wing systems reveals performance metrics for resonant flapping-wing flight

Flapping-wing insects, birds and robots are thought to offset the high power cost of oscillatory wing motion by using elastic elements for energy storage and return. Insects possess highly resilient elastic regions in their flight anatomy that may enable high dynamic efficiency. However, recent experiments highlight losses due to damping in the insect thorax that could reduce the benefit of those elastic elements. We performed experiments on, and simulations of, a dynamically scaled robophysical flapping model with an elastic element and biologically relevant structural damping to elucidate the roles of body mechanics, aerodynamics and actuation in spring-wing energetics. We measured oscillatory flapping-wing dynamics and energetics subject to a range of actuation parameters, system inertia and spring elasticity. To generalize these results, we derive the non-dimensional spring-wing equation of motion and present variables that describe the resonance properties of flapping systems: N, a measure of the relative influence of inertia and aerodynamics, and K^, the reduced stiffness. We show that internal damping scales with N, revealing that dynamic efficiency monotonically decreases with increasing N. Based on these results, we introduce a general framework for understanding the roles of internal damping, aerodynamic and inertial forces, and elastic structures within all spring-wing systems.


Introduction
Flapping-wing flight is one of the most energetically demanding modes of locomotion in nature and in engineered flying robotic systems. Actuators must provide power to overcome aerodynamic forces on the wings, generate inertial forces for oscillatory acceleration and deceleration, and counteract internal energy losses from imperfect power transmission [1]. If an oscillating wing is coupled to an elastic element such as a spring, the kinetic energy from the wing could be stored as elastic energy at the end of the wing-stroke and returned after stroke reversal. Many insects [1][2][3][4][5], birds [5,6] and even bats [7] have spring-like elements in the form of elastic materials in their thoraxes, muscles and tendons that may aid in reducing the energetic demands of flapping flight and improving flight efficiency (resonance) (figure 1a). However, the evidence that insects and birds operate near resonance largely relies on corelational observations of wingbeat frequency and wing inertia [9,10], or energetics arguments comparing metabolic and aerodynamic power [5,11,12].
If animals do rely on elastic energy storage for improved efficiency, then there are implications for the dynamics and energetics of those systems. One major example is that, to benefit from a spring, flapping wings must be actuated at a specific resonance frequency governed by the spring stiffness, body morphology and other factors such as aerodynamics and damping. Flapping at a higher or lower frequency leads to inevitable reduction in flight performance.
Early experiments on the relationship between wingbeat frequency and wing inertia provided compelling evidence that insects do oscillate their wings at resonance [9,10]. However, these experiments relied only on the manipulation of wing inertia (without accompanying measurements of thorax stiffness) and thus do not provide direct comparison of the spring-wing system's resonant frequency and the wingbeat frequency. Recently, there has been some effort to measure insect thorax elasticity and frequency response [13], but experimental limitations leave room for questions about whether insects are, in fact, flapping at a resonant frequency.
The impact of operation on or off resonance is directly related to how much a particular flapping-wing system benefits from the inclusion of a spring. A flyer with small wing inertia would have less excess kinetic energy to store and return in a spring than one with large wing inertia, but also would be less impacted by changes in wingbeat frequency. The question of where insects fall on this spectrum was first addressed by Weis-Fogh in his analysis of flapping flight efficiency [12]. He introduced the dynamic efficiency (η) as the ratio of aerodynamic work to the combined inertial and aerodynamic work, which serves as a measure of how much energy is expended on useful aerodynamic work versus wasteful inertial work. Weis-Fogh provided an analysis of the scaling of η in the absence of elasticity by introducing the non-dimensional variable N, which is the ratio of peak inertial force to peak aerodynamic force over a wing-stroke. In flapping-wing flight without springs, η was shown to monotonically decrease with increasing N, thus requiring larger energy input at larger N to sustain flight.
In subsequent text, we refer to N as the Weis-Fogh number to reflect his contributions to flight energetics.
Examination of the sinusoidal motion of a wing in the absence of elasticity reveals two sources of reaction force: inertial and aerodynamic forces. During a half-stroke of a wing, the inertial force associated with wing acceleration is at a maximum when the wing position is at its point of reversal, and inertial force decreases linearly with wing position and reaches zero at mid-stroke (figure 1b). At reversal, the wing speed is zero and thus the aerodynamic force at this point is zero, while at mid-stroke the aerodynamic force is maximum. Plotting these forces as a function of wing position (figure 1b) and normalized by peak aerodynamic force reveals that the inertial force has a maximum value of N, the Weis-Fogh number. Furthermore, the integration of these forces over the wing displacement provides the total inertial and aerodynamic work. As can be seen in figure 1b an ideal spring exactly matched to the inertial force of the wing would exactly cancel out the inertial work over a cycle in a parallel spring-wing system. In such a case the dynamic efficiency of the system would be η = 1 and the system would be operating at resonance. However, it is less clear how internal damping, frequency modulation, and different spring arrangements modulate the dynamic efficiency.
The primary mechanism of elastic energy storage in insects is resilin, a highly resilient elastic material first identified in patches of the locust exoskeleton, discovered and characterized in the 1960s by Weis-Fogh [3][4][5] (and subsequently identified in many other arthropods). It was shown to return greater than 97% of stored elastic energy, suggesting that insects have resilient components within their thorax that can facilitate efficient energy exchange and return. Thus, the historical choice not to include internal losses in the computation of dynamic efficiency appears to be a simplification based on the assumption that the losses due to aerodynamic forces are significantly larger.
However, recent experiments to characterize the energy storage and return in the hawkmoth (Manduca sexta) thorax have also demonstrated that small but significant energy loss occurs from internal structural damping [1]. Similar structural damping was also observed in cyclic oscillation of cockroach legs [14], possibly suggesting a general character of energy loss in exoskeleton deformations. Structural damping is a form of energy loss different from the more familiar velocity-dependent viscous damping. Materials that are structurally damped exhibit energy loss that is frequency independent and is instead governed by oscillation amplitude and elastic coefficient [15]. This is consistent with the interpretation of structural damping as a result of internal, microscopic friction that is not dependent on velocity. While the presence of highly elastic resilin suggests royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 18: 20200888 significant potential benefits from elastic materials, internal damping may preclude the energetic benefits of elasticity in spring-wing systems. Despite the more than 60 years of focus on resonance in insect flight, previous efforts at modelling or measuring spring-wing resonance in insects have fallen short by: (i) assuming that quasi-steady assumptions on aerodynamic forces in spring-wing resonance are valid, (ii) not including the effect of energetic losses from imperfect elasticity and (iii) focusing predominantly on the contributions of parallel system elasticity while disregarding contributions or limitations of series-elastic elements. Most importantly, we lack a common modelling and analysis framework to compare and contrast the energetic benefits of resonance across insects and man-made systems such as robots. Inspired by the original calculations of Weis-Fogh, we seek in this paper to develop a set of equations that govern the dynamics of parallel and series spring-wing systems using non-dimensional parameters that allow for comparative examination.
In experiments, we will measure how unsteady aerodynamic effects, specifically added mass and wing-wake interaction, influence the resonant behaviour of a flapping wing at hover. In order to achieve this, we compare simulations of a spring-wing oscillator subject to quasi-steady aerodynamic forces to a robophysical model of a flapping wing with known mechanical parameters subject to real fluid forces. Here we draw upon the work of van den Berg & Ellington [16], Sane & Dickinson [17], and others to use a flapping-wing robot that is dynamically scaled to that of flying insects. Unlike those studies, we do not prescribe the wing kinematics, instead using an elastic element in series between the wing and motor to observe resonant dynamics and emergent wing kinematics produced by varying actuation parameters. From our experiments, we will develop a model to understand how structural damping influences spring-wing dynamic efficiency using non-dimensional parameters. These efforts will provide a general understanding of how springs, wings, and body mechanics converge to enable energy efficient flapping motion as a function of morphology and wing kinematics.

Theoretical preliminaries: assumptions and motivation
To contextualize our study of spring-wing dynamics we first seek to outline the basic concepts of spring-wing systems. We will derive the equations of motion in the presence of aerodynamic and internal damping forces. We conclude this section with a non-dimensional representation of the equation of motion which produces two important parameters in spring-wing systems.

Undamped parallel and series wing-spring systems
The anatomies of flapping-wing animals feature a wide range of elastic element configurations that contribute to their flapping-wing dynamics. These arrangements can be expressed as combinations of springs in series and parallel with a moving inertial element, the wing (figure 1a). In both cases, the wing interacts with a time-and history-dependent nonlinear force from the surrounding fluid. To simplify the modelling of spring-wing systems, we consider the system as a one degree of freedom rotational joint (to emulate the wing hinge), where the joint angle θ is the wing angle along the stroke plane. If we neglect internal damping, the equation of motion for a parallel spring-wing system about this joint is where I t is the total system inertia, k p is the stiffness of the parallel elastic element, Q aero is the aerodynamic drag torque and T m (t) is the time-dependent torque applied to the wing. In the parallel configuration, the spring undergoes the same displacement as the wing and the muscle (actuator) acts directly on both the mass and spring. Nearly all spring-wing modelling of the insect thorax [12] and micro-aerial vehicles [18][19][20][21][22] has considered the parallel spring arrangement where muscle (actuator) is in parallel arrangement with the spring. In a series-elastic spring-wing system, a spring is placed between the actuator and the wing. Series elastic elements are well studied in vertebrate biomechanics as muscletendon units where a tendon is placed between the muscle and output [23]. Series elastic elements in flapping-wing flight may similarly be found in the muscle tendon units of birds [6,24]. In insects, series elasticity can arise from elastic tendons [25], elasticity in the wing hinge [4] or within the flight power muscle [26]. For simplicity of experiment design and to examine the differences between series and parallel systems, we analyse the series elastic spring-wing configuration. The equation of motion for a simple series elastic system may be written where the force acting on the wing arises from the displacement of input angle ϕ(t) relative to the wing angle θ (figure 2a). The difference between the angles is the deflection of the series spring with stiffness k s . When the system is at steady state (hovering), the series and parallel cases can be treated equivalently by rearranging equation (2.2) to reflect the parallel configuration:

Aerodynamic drag torque and added mass inertia
The wing experiences an aerodynamic resisting torque, Q aero , that opposes wing movement through the fluid. To make analysis of this system tractable, we will use a quasi-steady blade element estimate of aerodynamic torque consistent with previous quasi-steady methods for spring-wing systems [12] and micro-aerial vehicles [18][19][20][21][22]. Following the standard conventions for the quasi-steady blade element method, we express the aerodynamic drag torque The variable Γ is the drag torque coefficient and is a function of wing geometry (span, R; aspect ratio, AR; nondimensional radius of the third moment of wing area,r 3 ); pitch angle, α; drag coefficient, C D (α) and fluid density, ρ. royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 18: 20200888 In addition to aerodynamic torque, the acceleration of a wing within a fluid leads to an additional inertia; an 'added' or 'virtual' mass, as it is sometimes called [27]. We use a method from [27] for modelling the mean added mass inertia, I A (see electronic supplementary material, §11.1). The total inertia is computed as I t = I sys + I A , where I sys is the inertia of the wing and wing transmission. In the insect flight system, the wing hinge acts as a mechanical transmission, converting linear muscle actuation to angular wing motion. It is possible that additional inertial terms from the reflected inertia of the oscillating muscle and thorax may be important. However, there is little known about the specific motion and inertia of the wing hinge transmission and the role of muscle and thoracic inertia, so we will disregard their effects in this paper.

Structural damping in the insect thorax
Recent experiments to measure the damping response of the hawkmoth (Manduca sexta) thorax [1] and cockroach (Blaberus discoidalis) leg joints [14] have both identified structural or hysteretic damping as the dominant source of energy loss. Consistent with these observations, we seek to consider the effects of structural damping on spring-wing system dynamics. Structural damping is a common source of energy loss in biomaterials [28] that differs from viscous damping in that there is no velocity dependence in the structural damping force. For general oscillatory motion, the structural damping force can be included as a modification to the spring constant, K = k(1 + iγ), where k represents either the parallel or series spring. The coefficient γ is the structural damping loss factor [15] which has been found to be γ = 0.2 for cockroach leg joints [14] and γ = 0.1 for the hawkmoth thorax [1]. For constant sinusoidal motion at a single frequency, ω, the structural damping force can be represented as a viscous-like force with a coefficient that scales with frequency where the angular velocity _ b is the relative speed of spring compression (for parallel systems . The presence of ω in the denominator makes the structural damping force frequency-independent, unlike typical viscous damping (see electronic supplementary material, §11.2, for derivation).

An organizational framework for spring-wing systems
In this final section of the motivation, we introduce a set of non-dimensional variables that govern general springwing dynamics. As discussed in the Introduction, the Weis-Fogh number N is the ratio of maximum inertial force compared to maximum aerodynamic force over a cycle. Assuming sinusoidal wing motion with amplitude θ 0 and frequency ω, the maximum inertial torque is I t θ 0 ω 2 and maximum aerodynamic torque is Gu 2 0 v 2 resulting in the Weis-Fogh number This quantity should dictate how important a role spring elements can play in energetic efficiency at hover. For N < 1, aerodynamic forces dominate and kinetic energy may be fully dissipated into the surrounding fluid over each wing stroke; no elastic storage is needed. However, for N > 1, excess kinetic energy from the wing can be recovered by a spring. Observations from biological and robotic flappingwing flyers indicate that N roughly varies between 1 and 10 for a broad range of insects [12]. In order to compare across insect species, we non-dimensionlize the dynamics, assuming the wing oscillates sinusoidally at a frequency ω, and with amplitude θ 0 . We begin with the full dynamics equation for the parallel spring-wing system including structural damping: We introduce dimensionless wing angle q and dimensionless time τ and substitute into equation (2.7) (see electronic supplementary material, §11.3). We arrive at the following dimensionless equation of motion for the spring-wing system: where we have defined the reduced parallel stiffness, and the Weis-Fogh number, N, is in the coefficient of the aerodynamic torque. The normalized torque in the parallel system isT Performing a similar substitution for the series system we arrive at the equation below: where the normalized torque in the series system iŝ We provide a full derivation of these equation in the electronic supplementary material, §11.3.
Through the change of variables in the parallel and series spring arrangements we have arrived at two nearly identical non-dimensional dynamics equations in equations (2.8) and (2.11). The differences between the series and parallel systems in this form are all contained in the actuation variables,T s andT p . Thus, while the forced actuation of these systems may result in different dynamics, the similar structure of the non-dimensional dynamics equations indicates that in both systems three variables likely govern the behaviour: N, K and γ. As an exemplary demonstration of this, if either the series or parallel system is driven to steady-state oscillation and then the external actuation is turned off (T = 0 or ϕ = 0), both equations (2.8) and (2.11) become equivalent. The objective of this paper is to seek to understand how N, K and γ influence dynamic efficiency and resonance of a simple series spring-wing system.

Robophysical experiment design
We developed an experimental spring-wing system to study flapping-wing behaviour in the context of realistic fluid forces (figure 2). The quasi-steady modelling presented in the previous section greatly simplifies the unsteady, history-dependent aerodynamic phenomena involved in flapping flight. A dynamically scaled physical wing serves as a reference against which we can evaluate the performance of the quasi-steady model. The system consists of a servo motor capable of accurate position control, a moulded silicone torsional spring with linear elasticity and structural damping, and a simple fixed-pitch wing element attached to a rotary shaft and submerged in a 115-gallon plastic tank filled with water (30" × 30" × 30", Chem-Tainer). The tank was selected to be large enough to minimize fluid interactions with the side walls, floor and water surface. The wing is situated near the centre of the tank, such that the wing is always at least 10 wing-chord-lengths from the walls and floor and 5 wing-chord-lengths from the surface of the water. This is consistent with other studies of flapping wings that use water as a working fluid [29][30][31]. See electronic supplementary material, §11.4, for photos of components and further details.

Motor selection and system friction reduction
A high-torque servo motor (Teknic Clearpath SDSK) was chosen to drive the system under closed loop angular position control. The servo is able to provide substantially more torque than that experienced by the wing in the fluid, effectively decoupling the motor and wing dynamics. We monitor the motor and wing angle using two optical encoders (US Digital, 4096 CPR). To reduce the influence of friction on the wing motion, we used two radial air bearings (New Way, S301201) which resulted in negligible bearing friction. The shaft was supported vertically by an axial thrust bearing, which did contribute a small amount of friction.

Reynolds number scaling
To ensure that we match the aerodynamic behaviour of small insect wings we chose experimental parameters to dynamically scale our system. Consistent with previous dynamically scaled experiments [17,32], we sought to maintain a Reynolds number in the range of that of small flapping-wing insects, Re = 100-10 000. We define Reynolds number based on standard methods, using wingtip velocity as the flow speed and wing chord as the characteristic length [32]. We choose water as a working fluid (ρ = 997 kg m −3 ) and choose wing geometry (rectangular, 10 × 3.6 × 0.5 cm) and a range of actuation parameters (10-64 deg amplitude, 0.5-4.1 Hz frequency) where the resulting wing kinematics had Re ≈ 10 3 − 10 4 . Note that since the wing amplitude is an emergent property of the system due to the series spring configuration, so too is the Reynolds number of an individual test.

Weis-Fogh number scaling
In order to achieve a wide range of Weis-Fogh number (equation (2.6)) in experiment, we must be able to adjust the ratio of peak inertial torque to peak aerodynamic torque. We chose to keep the drag torque coefficient Γ constant, so we vary N by changing wing inertia and wing amplitude θ 0 . Wing inertia I t was varied by adding acrylic and aluminium discs to the wing shaft, leading to total inertias of I t = [10.5, 15.9, 23.2] × 10 −4 kg m 2 . Our experiments resulted in a range of N between 1 and 5 when operating at the resonant frequency, consistent with many insects [12].

Aerodynamic calculations
We used a rigid rectangular wing with a fixed vertical pitch (α = 0), which simplifies modelling and motor control, and eliminates any energy storage and return from a flexible wing. The drag torque coefficient for a rectangular wing from equation (2.4) is Γ = 1.07 × 10 −3 kg m 2 , where the coefficient of drag is constant at C D (0) = 3.4 (taking this value from [17]). We also compute added mass inertia, I A , for a calculated value of 3.465 × 10 −4 kg m 2 (see electronic supplementary material, §11.1, for derivation).

Silicone spring fabrication and evaluation
We used custom fabricated silicone torsion springs for the series spring element (figure 3a). Silicone was chosen because it can be cast into custom shapes and has a linear elastic response over large strain [33]. The springs were designed royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 18: 20200888 with a cylindrical profile with flanges on each end to facilitate coupling to the motor and wing shafts (figure 3a). Detailed information about the design and fabrication process may be found in electronic supplementary material, §11.6.
We used three spring designs with torsional stiffness values of K s ¼ [0:163, 0:416, 0:632] N m rad À1 . Figure 3 shows the results of experiments to characterize the spring mechanical properties. We subjected springs to both cyclic and static loading conditions (figure 3a,b) to measure their elastic and damping properties. The spring torque response was linear over the range of angles tested (figure 3b,c) with stiffness values that are consistent with the predicted torsional stiffness (figure 3d ). In dynamic testing, we observed a small amount of hysteresis in cyclic loading experiments indicating the presence of damping within the spring (figure 3b). There is evidence that silicone rubber has a combination of viscous and structural damping [34,35], so we set out to test whether a viscous or structural model fits best. Following a similar procedure as [1], we oscillated each spring sinusoidally across frequencies between 1 Hz and 10 Hz and amplitudes of 10, 20 and 30 degrees. We measured the angle-torque relationship during these tests and fitted a viscous and a structural damping model to the data for each spring. We expected that the model that fitted best would have the least variation in the fit coefficients across the range of frequencies. We computed the mean and standard deviation of the fit coefficients across the range of frequencies and found the standard deviation of the coefficients as a percentage of the mean coefficient. The standard deviations in the viscous model coefficients were 74.8%, 83.2% and 62.1% of the mean for springs K1, K2 and K3, respectively, whereas those of the structural model were only 11.7%, 6% and 22.3%, respectively. Therefore, we chose to use the structural model to model the damping, with the knowledge that the true damping is more complex but outside the scope of this paper. We measured the structural loss moduli shown in figure 3e. For comparison, the loss modulus is γ = 0.1 in hawkmoths [1] and γ = 0.2 in the legs of cockroaches [14]. Each test was run for 30 s, during which encoder position readings were recorded from both the motor and wing by a NI PCIe board (6323 Multifunction IO Device), sampled at 1 kHz. Encoder readings were saved as text files in MATLAB and processed.

Simulation
We developed numerical simulations that modelled the parallel and series configurations of the spring-wing system based on equations (2.1) and (2.2). The parameters of the simulations were based on measured and calculated parameters from the experiment. Simulations were performed in MATLAB using the ODE45 numerical integration function. The integration algorithm performs an adaptive step integration with absolute and relative error tolerances set at 10 −3 and 10 −6 , respectively.

Data analysis
We performed identical analysis of both simulation and experimental data. We collected wing and motor angle data, computed wing and motor angular velocity, and found the amplitude, frequency and relative phase of each by fitting the data to functions θ(t) and ϕ(t), respectively: See electronic supplementary material, §11.7, for details. In order to identify the resonant peak, we define a nondimensional term, kinematic gain, which is the ratio of motor angle amplitude and wing angle amplitude: For a resonant system, kinematic gain is maximum at the resonant frequency, ω r ( figure 4a,b). We also compute the quality factor of the oscillator dynamics, Q, using the following definition:

4)
where ω r is the resonant frequency and Δω is full width at half maximum. The quality factor is a metric of the sharpness of the resonant peak: high Q means kinematic gain G K drops royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 18: 20200888 off as ω moves away from the resonant frequency, while low Q means G K changes slowly with varying ω.
In order to identify ω r for a set input amplitude experiment, we locate the frequency that maximizes kinematic gain. We fit a fifth-order polynomial to the 12 points closest to the peak measured G K to get a smooth approximation of the resonance curve. The maximum value of the polynomial is the peak gain, and the frequency corresponding to it is the resonant frequency. When reporting N, dynamic efficiency, etc. 'on resonance', we use the experimental configuration with a frequency closest to the resonant frequency. As a result, some nominally resonant points are not exactly on the resonant peak, but may be off by as much as 0.1 Hz.
We use the position measurements, their derivatives, and the known mechanical parameters of the system to estimate the torques on the system: aerodynamic, inertial, elastic, motor and structural damping. and Note that we compute the equivalent non-dimensional torques and kinematics using the terms defined in equations (2.8) and (2.11). Lastly, we compute the dynamic efficiency of the oscillatory motion. Weis-Fogh's definition of dynamic efficiency from the introduction neglects elasticity and internal damping, but for our purposes we need dynamic efficiency to take into account aerodynamic, inertial, elastic and damping work. A simple way of doing so is to define the dynamic efficiency as the ratio of aerodynamic work to total work done by the actuator: Note that in the motor work integral R(x) is a rectification function defined as R(x) = x for x > 0 and R(x) = 0 for x < 0. This rectification function accounts for the fact that negative mechanical work in insect flight muscle does not contribute significantly to metabolic energy expenditure [12,36].

Kinematic gain varies with actuation and system properties
We performed 3249 tests with varying input parameters (amplitude and frequency) and mechanical system parameters (spring stiffness and inertia) and measured the emergent wing kinematics (figure 5). Results for the nine inertia and stiffness combinations are shown as heatmaps with colour indicating kinematic gain. The arrows on the right and top of the figure denote the directions of increasing stiffness and inertia, respectively. For all stiffness and inertia combinations, we observed that the peak G K occurred at the lowest actuation amplitudes. At low amplitudes, aerodynamic damping is minimized and less energy is lost to the surrounding fluid, allowing for higher kinematic gain. The maximum kinematic gain increased with increasing system inertia in all cases, reaching a maximum when the system is a combination of a soft spring and high inertia (bottom right corner).
The resonant frequency decreased as the motor amplitude (and thus the emergent wing amplitude) increased. The dashed purple line in figure 5 shows peak kinematic gain for each motor amplitude. This decrease in resonant frequency is consistent with simulation predictions shown as solid black lines in figure 5. We discuss this model in the first discussion section below.

Flapping resonance with quasi-steady aerodynamics
Following the experiments, we sought to see how much of the observed dynamics is predicted by the simplified, quasi steady equations of motion described in §2. The real aerodynamic loads on the wing are time-and history-dependent, so it is not clear to what degree those unsteady loads affect the resonance properties of the system at steady state.

Natural frequency of the system matches lightly damped resonant frequency
At low motor amplitude, the system dynamics are only weakly affected by the aerodynamic force and thus the resonant frequency should be determined by the natural frequency of the spring-wing system. Comparison of the experimentally measured resonance frequency at the lowest amplitude with natural frequency computed from spring stiffness and inertia displays extremely good agreement (figure 6a). Thus our intuition is confirmed that low amplitude wing oscillation leads to small aerodynamic damping and the system better resembles a weakly damped spring-mass oscillator.

Quasi-steady simulation predicts experimental resonant frequency
As motor and wing amplitude increase, we observe reductions in both kinematic gain and resonant frequency. This behaviour is consistent with the quadratic velocity damping generated by the aerodynamic force on the wing and has been observed in other spring-wing experiments [20]. We sought to determine if a quasi-steady aerodynamics model was sufficient to model the spring-wing behaviour observed in experiment. By performing numerical simulations with identical parameters as the experiments of figure 5 we measured the simulated resonant frequency across all amplitudes. We find that the quasi-steady simulation resonant frequency agrees well with the experimental resonance frequency ( figure 6b). This suggests that quasi-steady aerodynamic modelling is sufficient to capture the spring-wing relationship between amplitude and resonance frequency across all motor amplitudes.

Quasi-steady simulation over-predicts kinematic gain
Despite the agreement in resonant frequency, we observed a difference between the kinematic gain in the experiment and the simulation (figure 6c). The simulation over-predicted the kinematic gain for all resonance experiments by a maximum of 20% at the highest experimental gains. The over-prediction error grew with increasing gain, suggesting a systematic error in the modelling of the spring-wing system. The combination of disagreement between simulation and experiment in kinematic gain and the good agreement in resonant frequency suggests to us that unmodelled dissipation from system friction is likely the cause. Coulomb friction in the bearings would not influence the system resonant frequency, but would decrease the output amplitude, consistent with our observations. We opted not to include friction in our modelling for two reasons: (i) we kept only the biologically relevant damping terms in the system equations and (ii) modelling friction can be complicated due in part to highly nonlinear stick-slip phenomena [37].

Dynamic efficiency is amplitude and frequency dependent
To determine the efficacy of elastic springs for energy reduction in a flapping-wing system, we calculated the dynamic efficiency, η, across all system and actuation parameters ( figure 7). Generally, we observe that dynamic efficiency is maximum along the resonance curves for all experiments. These results are consistent with the interpretation that maximum kinematic gain corresponds to maximum energetic benefits of having a spring, e.g. that the system is at resonance. Notably, the dynamic efficiency is very sensitive to oscillation frequency for low motor amplitudes while higher motor amplitudes show a very broad dynamic efficiency. The results at high amplitude are consistent with the broad dynamic efficiency versus frequency curves measured in experiments on a flapping-wing robot in [19].

Discussion
In this discussion, we recall the preliminary framework established in §2. Through a change of variables we were able to express the equations of motion of the series and parallel spring-wing systems (equations (2.1) and (2.2)) with nearly identical expressions. In this discussion, we now seek to interpret the series elastic experiment and simulation results above in the context of the non-dimensional variables: the reduced stiffnessK, the Weis-Fogh number N and the structural damping coefficient γ.

The Weis-Fogh number N governs resonant properties in spring-wing systems
The frequency response and kinematic gain of the springwing system indicate that the resonant behaviour of the system is dependent on the flapping amplitude and frequency. In a standard spring-mass system with a viscous damper, frequency-dependent kinematic gain is to be expected. However, the dependence of the system resonant royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 18: 20200888 frequency on oscillation amplitude in the series spring-wing system (figure 5) is different from that of the standard viscously damped spring-mass system (unless the damping is close to a critical value). Similar to the way the springmass system is often reduced to non-dimensional ratios (damping ratio ζ and reduced frequency ω/ω n ) that govern the oscillation behaviour, we here seek to show how the non-dimensional parameters of equations (2.8) and (2.11) govern the resonance properties of the system.

N is a measure of aerodynamic damping in spring-wing resonance
An important metric of a resonant system is the quality factor, Q, a measure of how damped the oscillator is and which, in experiment, is the sensitivity of kinematic gain to frequency change. For a linear spring-mass system, Q is independent of oscillation amplitude, but it is inversely proportional to the damping ratio, ζ. In the non-dimensional spring-wing equations, the two terms that govern system damping are the Weis-Fogh number N (the prefactor of the aerodynamic term) and the structural damping loss modulus γ. We expect that for any useful spring-wing system, the energy loss due to aerodynamic damping should be much larger than the parasitic energy loss from internal structural damping. Thus, we are inspired to examine how quality factor varies with the Weis-Fogh number.
Examining the dependence of Q on the Weis-Fogh number, we observe that the quality factor grows approximately linearly with the Weis-Fogh number (figure 8a) across the ranges observed in experiment. Simulations further reveal a linear relationship between Weis-Fogh number and system quality factor given by where c 1 = 0.18 ± 0.01 and c 2 = 1.11 ± 0.03 (95% CI). By analogy with the linear spring-mass system, the linear relationship between Q and N reveals that the quantity 1/N is comparable to the damping ratio for linear spring-mass systems. The importance of this relationship is that Q has historically been used to evaluate the ability of insects to benefit from elastic energy storage. Ellington, using a version of Q defined as the ratio of peak kinetic energy versus energy dissipated per stroke, has reported Q values for a variety of insects such as the fruit-fly Drosophila melanogaster, hawkmoth Manduca sexta, and the bumblebee Bombus terrestris [38]. Ellington uses those values to argue that insects would benefit from resonant dynamics. While it is not possible to directly compare the values he reports due to the use of different expressions for Q, we can use the functional relationship in equation (5.1) to also estimate quality factors for the fruit-fly (N ≤ 1, Q ≈ 1.2), the hawkmoth (N = 3.6, Q ≈ 1.7), and the bumblebee (N = 3.1, Q ≈ 1.6). Weis-Fogh numbers for the relevant animals were provided by Weis-Fogh in [12]. Quality factors above 1, in this case, confirm

The resonant frequency of a series spring-wing system varies with N
Examination of the Weis-Fogh number shows that N grows inversely with wing amplitude (equation (2.6)). Thus, experiments that were performed at low motor amplitude correspond to spring-wing systems with large N. This observation, coupled with the insight from above, immediately provides understanding for why the resonant frequencies at low amplitudes match the system's natural frequency (figure 6a). At low amplitude (high N), the influence of damping is very small (scales as 1 N ) and thus increasing N in spring-wing systems results in minimizing the effects of aerodynamic nonlinearity and energy loss.
However, this insight does not fully explain why the resonance frequency should change with amplitude. Examination of the relationship between the non-dimensional variables that relate frequency (K) and damping (N) at resonance shows a tight locus of points that indicate a potential relationship (figure 8b). To understand this relationship, we follow the method originally introduced by Bennett et al. for the analysis of series elasticity in whale flukes [39]. In that paper, the authors determine the instantaneous aerodynamic, inertial, and elastic power contributions in the system and then minimize the total power consumption over a wing-stroke (see electronic supplementary material, §11.8) to identify an optimal combination of parameters. The combination of parameters that minimize total power is given by the following relationship: We plot this relationship in figure 8b and observe good agreement (R 2 = 0.868, RMSE = 0.22) between the predicted optimal series spring-wing relationship and observed relationship at resonance. The observed differences between experiment and theory could result from measurement error, or alternatively friction, structural damping, or aerodynamic effects that are unmodelled in the original derivation.
However, solving equation (5.2) in terms of the actual resonant frequency ω r and system parameters (see electronic supplementary material, §11.8) shows very good agreement over the frequency, motor amplitude parameter space (figure 5, black line). For completeness, we note that for a parallel spring-wing system, the resonance relationship is constant:K ¼ 1 for all N [39]. The overall intuition from this resonance analysis is that the Weis-Fogh number plays an important role in determining the spring-wing system's resonant behaviour. Furthermore, the relationship between N andK at resonance implies that these two variables likely define a non-dimensional parameter space for general spring-wing systems of arbitrary morphology, mechanical properties and actuation. In the following sections, we will expand that space by exploring how N,K and structural loss factor γ influence dynamic efficiency.

Energy loss from structural damping scales with
Weis-Fogh number in real spring-wing systems Recent discoveries that structural damping in hawkmoth thoraxes [1] and cockroach legs [14] may cause significant energy loss prompts us to investigate how structural damping contributes to spring-wing energetics. To gain insight into energy loss scaling, we consider the parallel arrangement of a spring and a wing ( figure 1a). The parallel arrangement, as opposed to the series arrangement, is convenient because the spring (and structural damper) is subject to identical displacements as the wing hinges, and thus all forces (spring, damping, aerodynamic and inertial) can be represented as functions of wing angle. In the series case, the actuation and wing trajectories are out of phase and therefore more complicated, if not impossible, to express analytically. In the parallel system, we can conveniently visualize all the relevant forces acting on the wing by plotting torque contributions versus wing angle (figures 1a and 9a). We now consider the scaling of each of these forces with the non-dimensional system parameters. Following the method introduced by Weis-Fogh, we normalize all the non-dimensional torques (equation (2.8)) by the peak aerodynamic torque at mid-stroke. The result is the following set of non-dimensional torques, expressed as functions of royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 18: 20200888 normalized wing angle, q. The tilde symbol represents torques normalized by peak aerodynamic torque: See electronic supplementary material, §11.9, for the full derivation of these terms and §11.10 for an additional derivation of the torque due to viscous damping. In the parallel system, the ideal spring torque is exactly opposite the inertial torque and thus cancels the inertial work throughout the wing stroke,Q elastic ¼Q inertial . This relationship implies thatK ¼ 1 in the ideal parallel system. Converting this expression back to dimensional form returns the expected relationship that defines a parallel resonant system, k p = I t ω 2 .
Critically for analysis of spring-wing energetics, all of these equations can be integrated over the wing stroke to provide expressions for the non-dimensional work the wing has performed. In the case of parallel resonance (K ¼ 1), we simply need to integrate theQ aero andQ structural terms since the inertial and elastic work exactly cancel. Performing these integrals over the range q = [ − 1, 1] results in the following expressions: Since aerodynamic and structural damping work are the only sources of energy loss (and thus the only required energy input) in the system we can now express the dynamic efficiency in terms of these two energy losses: h ¼W aerõ W aero þW structural (5:6) We now have an expression for the dynamic efficiency of a parallel spring-wing system operating at the resonance frequency. The dynamic efficiency is a function of only the Weis-Fogh number, N, and the structural damping γ. Examination of this expression indicates that if there is no structural damping in the system (γ = 0) the dynamic efficiency is constant, η = 1, across all N, indicating that all input work goes towards aerodynamic work at steady state. However, for any non-zero structural damping in the transmission, the dynamic efficiency is a monotonically decreasing function of N, since a portion of the input work must be diverted to overcome the internal structural damping.
This analysis of a parallel spring-wing system has provided insight into how structural damping influences the mechanical efficiency of the flight transmission, i.e. dynamic efficiency. We would like to consider the same theoretical analysis for a series spring-wing system, but in this case the theoretical approach becomes intractable. In a series springwing system, the relative phase difference between wing and actuator changes with actuation parameters and thus the relative spring extension (and rate of extension) is not as easily determined. Bennett et al. [39] presented an elegant reformulation of the problem to generate a closed form expression for the actuator power required in a series spring-wing system. However, this method cannot be used to include structural damping because it would require knowledge of the actuator input kinematics. Thus, to determine the influence of structural damping in a series springwing system, we will resort to numerical methods in the next section.

Structural damping reduces peak dynamic efficiency in series spring-wing systems
The dynamic efficiency η for series spring-wing systems was determined through numerical simulation of the non-dimensional equation of motion (2.11) across a range of N andK. As expected in the case of no structural damping, γ = 0, the dynamic efficiency was η = 1 for all values of N andK along the series resonance curve (equation (5.2)) (figure 10a). These simulations also allow us to observe how η varies across the full (N,K) parameter space. In general, we observe stronger sensitivity of η to changes inK as N increases. Recalling our analysis of the connection between quality factor and Weis-Fogh number N, it is clear that, for a spring-wing system with fixed morphology and no structural damping, as N increases, the variation of η with frequency becomes more significant.
For systems with non-zero structural damping, as expected in any real system (and as observed recently [1,14]), the dynamic efficiency in general decreases with increasing N for all values ofK. For small structural damping (γ = 0.1; figure 10b) we observe a general similarity of the dynamic efficiency to the undamped case. The gradient of the dynamic efficiency can be observed by the spread of the contour lines, where a steep gradient is indicated by closely spaced contour lines. In the undamped case, the resonance relationship follows the line of minimum gradient in dynamic efficiency (minrh): the resonance curve exactly follows the contour line of η = 1 (and thus rh ¼ 0). Examination of η for increasing structural damping indicates the curve of minimum rh likely differs from the undamped resonance as illustrated in figure 10a, where the dashed lines represent estimates of the line of minimum gradient to guide the eye.
Evaluating the dynamic efficiency across the different values of γ shows that η monotonically decreases with N, consistent with our analysis of the parallel spring-wing system. In the inset of figure 10b, we show the dynamic efficiency results at resonance for the experimental series spring-wing system. The experiment exhibits the same monotonically decreasing trend of η with N. However the magnitude of η in the experiment differs from that predicted by the simulation with structural damping alone. Thus, similar to the previous comparison of simulation and experiment (figure 6) we observe qualitatively similar trends between experiment and simulation; however, the experiment exhibited lower η, likely due to additional sources of energy loss from friction and other unmodelled effects.
Overall, both the experimental and simulation results provide evidence that, for any spring-wing system with structural damping, the dynamic efficiency decreases monotonically with increasing N. Furthermore, dimensional analysis of the viscous damping model (see electronic supplementary material, § 11.10) suggests that internal damping of many types has the same effect on dynamic efficiency. These results are a bit counterintuitive from the discussion above, where quality factor is observed to scale linearly with N (figure 8a). For a perfect spring-wing system (γ = 0) it is true that increasing N diminishes the relative energy loss from aerodynamic work compared to the total energy of the oscillator (the definition of the quality factor). However, in the presence of internal energy losses, the actuators have to do extra work to overcome internal body damping, which scales with N. With internal damping, the quality factor might still be the same (since it is a ratio of energies) but the internal damping decreases the useful energy of the oscillator and thus dynamic efficiency.

Intermediate Weis-Fogh number may balance damping losses and elastic benefits
The scaling of dynamic efficiency with N may help explain why flapping-wing insects and robots tend to have an N in the range of 1-10 across three orders of magnitude in body mass (figure 11). A very low N means that aerodynamic forces dominate and wing kinetic energy will be dissipated by aerodynamic work rather than be stored in a spring [11]. High N, however, increases the internal energy losses from structural damping and reduces the benefits of elastic energy storage. Thus, flapping-wing insects and birds with Weis-Fogh number in an intermediate range may balance the positive benefits of spring-wing energy exchange with the parasitic energy losses of internal structural damping. In order to achieve high dynamic efficiency at hover, wing geometries, flapping amplitudes, and wingbeat frequencies may be tuned to maintain operation in this restricted regime of Ns.
In addition to the constraints of damping, high N biological spring-wing systems may not be possible due to biomaterial or physiological constraints. For example, high N spring-wing systems that are capable of hover (implying modest wingbeat amplitude, θ 0 ) would require large wing inertia which may be impractical due to the possible imposition of extra system weight. Similarly, from the elastic materials perspective a high-N spring-wing system with large wing inertia would require extremely stiff and resilient elastic elements to operate at resonance. It may be that for extremely large N spring-wings the required elastic stiffness may exceed the practical regime of biological materials. Both of these considerations however do not necessarily limit robotic systems from being developed with high N through appropriate inertial and stiffness design considerations.
Lastly, it is unlikely that dynamic efficiency at hover is the only factor that dictates this range of morphologies. Other factors, such as the effect of N on transient dynamics and control of wing kinematics, are likely to be significant. For  Figure 11. The Weis-Fogh number for flying insects and robots as a function of body mass (left plot) and in fixed flapping-wing robots (right boxplot). Biological data are reported from [12], and robot data from [19,[40][41][42]. Most observed insects and robots lie within the mass-independent range of N ∈ [1, 10].
royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 18: 20200888 example, in a spring-wing system with high N the wing dynamics are dominated by inertial effects and thus wing kinematics are likely insensitive to transient changes in aerodynamic forces. Such a high-N system may have wingbeat kinematics that are relatively stable in the presence of gusts of wind. However, in the opposite case of a low-N springwing system, the flyer would be able to more easily modulate wing kinematics and possibly wingbeat frequency for control purposes. In these cases and others, the Weis-Fogh number may provide a baseline for insight in broad comparative study of flying insects, enabling identification of commonalities between species as well as exceptional cases that merit further study.

Conclusion
Many flapping-wing insects and birds possess elastic elements in their body that may reduce the power demands of flapping-wing flight. However, recent experiments have demonstrated that insects are also subjected to internal power loss from the deformation of their thorax. In this paper, we have introduced three non-dimensional variables for general spring-wing systems that govern oscillatory behaviour and dynamic efficiency. Inspired by the foundational work of Weis-Fogh, we re-introduce the ratio of maximum inertial force to aerodynamic force as the Weis-Fogh number, N. Experiments and simulation illustrate that N is a fundamental parameter of spring-wing systems, analogous to the quality factor of a linear spring-mass system. However, when spring-wing systems have internal structural damping, we observe that dynamic efficiency decreases with increasing N on resonance, reducing the potential for useful energy storage and return. Overall, these results provide a generic framework to understand spring-wing systems which may enable us to learn more about the inter-relationships of morphology and actuation in flapping-wing insects and birds.