Abstract

Flapping-wing flight is a challenging system integration problem for designers due to tight coupling between propulsion and flexible wing subsystems with variable kinematics. High fidelity models that capture all the subsystem interactions are computationally expensive and too complex for design space exploration and optimization studies. A combination of simplified modeling and validation with experimental data offers a more tractable approach to system design and integration, which maintains acceptable accuracy. However, experimental data on flapping-wing aerial vehicles which are collected in a static laboratory test or a wind tunnel test are limited because of the rigid mounting of the vehicle, which alters the natural body response to flapping forces generated. In this study, a flapping-wing aerial vehicle is instrumented to provide in-flight data collection that is unhindered by rigid mounting strategies. The sensor suite includes measurements of attitude, heading, altitude, airspeed, position, wing angle, and voltage and current supplied to the drive motors. This in-flight data are used to setup a modified strip theory aerodynamic model with physically realistic flight conditions. A coupled model that predicts wing motions is then constructed by combining the aerodynamic model with a model of flexible wing twist dynamics and enforcing motor torque and speed bandwidth constraints. Finally, the results of experimental testing are compared to the coupled modeling framework to establish the effectiveness of the proposed approach for improving predictive accuracy by reducing errors in wing motion specification.

Introduction

Flapping-wing air vehicles (FWAVs) are a unique set of unmanned aerial vehicles (UAVs) that derive inspiration from natural sources like birds, bats, insects, and even marine life. Flying animals achieve a broad range of capabilities from hovering and precision maneuvers to long-range cruising and efficient soaring. Due to this potential for excellent versatility, flapping-wing flight may provide a useful alternative to airplanes and rotorcraft. Flapping-wing aerial vehicles may be grouped into two main size categories, corresponding to the sources of biological inspiration: insect and avian. Due to favorable scaling of aerodynamics with larger size, avian-scale flapping-wing air vehicles are capable of lifting greater loads for longer periods of time [1]. Hence, avian-scale flapping-wing air vehicles are currently better suited for missions involving external load carrying, longer endurance flight, and operation in adverse conditions. This study is focused on avian-scale flight, which may be characterized as forward flight, rather than hovering flight often exhibited by insects.

Typically, avian-scale flapping-wing aerial vehicles rely on the interaction of a few key components to support flight [2]. First, a propulsion system typically consisting of a motor and gear train provides the periodic wing plunging motion necessary to generate forces. Second, compliant wings are flapped and either actively or passively deformed in response to actuator inputs or aerodynamic and inertial loading, resulting in useful lift and thrust production. Third, a source of energy provides the ability to freely fly for a certain period of time. In addition to these key components, most vehicles consist of several additional parts including control systems and structural support. For designers, the interaction of these key subsystems presents a significant challenge, due to the interdependence of the performance of each subsystem on vehicle capabilities.

Historically, successfully flying avian-scale flapping-wing aerial vehicles have been designed using experimental performance characterization. A variety of flapping-wing aerial vehicles have been created with this approach and have clearly demonstrated the breadth of capabilities that may be realized [37]. By anchoring the vehicle to a load cell, researchers at the University of Maryland, College Park, MD have examined the scaling of forces with flapping rate in a small vehicle [8]. Hubel and Tropea studied the effects of varying Reynolds numbers in a wind tunnel experiment to explore the appropriateness of quasi-steady assumptions in modeling [9]. Shkarayev et al. [10,11] and Song et al. [12] investigated flexible membrane wings and discovered the sensitivity of lift and thrust production to wing membrane pretension. Each of these efforts realizes desired functionality through a direct experimentation approach but does not explicitly address the nature of wing and motor interactions.

A significant issue with experimental studies of force production under varying flapping conditions is the rigid mounting to a force transducer. This source of error in flapping-wing experiments was explored by Caetano et al. in the development of the Delfly flapping-wing aerial vehicle [13]. By restricting the motion of the flying vehicle for testing purposes, natural body motions that arise in response to flapping are prevented which causes a change in the measured forces [13]. A strategy to mitigate this effect is to equip a freely flying vehicle with instrumentation to avoid altering the natural body response to flapping. To date, there is only one example of in-flight testing of a flapping-wing aerial vehicle with an instrumentation suite, conducted by Grauer and Hubbard [14]. This test was able to collect some basic in-flight data, but suffered from noisy sensor data that required time-averaging of results and lacked measurements of power consumption, vehicle position, altitude, and airspeed.

Aerodynamic modeling provides an alternative to experimental development efforts. Aerodynamic models are used to generate performance predictions for flapping-wing aerial vehicles in a variety of flight conditions by specifying wing motions and body orientation, then solving for lift and thrust forces. DeLaurier developed a well-known model of flapping-wing flight based on modified strip theory, which assumes quasi-steady aerodynamics and sums forces at each location and time step to generate performance estimates [15]. DeLaurier's work relies on Theodorsen's plunging airfoil model to account for force production due to leading-edge suction forces [16]. The method is simple to setup and run, but showed limited accuracy when applied to wings that exhibit larger flapping amplitudes by Mazaheri and Ebrahimi [17]. Hunsaker and Phillips have compared Theodorsen's model to computational fluid dynamics simulation results, revealing that the basic Delaurier model formulation provides reasonable results in lower flapping frequencies [18]. A dynamic stall approximation by Kim et al. exhibits improved prediction accuracy for cases of larger angles of attack [19]. Smith et al. presented a compact review of aerodynamic methods for flapping wings, which includes a useful comparative discussion of the merits of several approaches. In addition, they present a comparison between a lower fidelity quasi-steady method and an unsteady vortex panel method to highlight the strengths of each strategy [20]. A common feature of the models reviewed here is that flapping motions are specified without consideration of the interactions between the drive motors and the predicted forces arising due to aerodynamic and structural loads. This assumption implies a motor with sufficient torque and power to achieve specified motions, which is often an invalid assumption in real operation.

In addition to modeling aerodynamics, researchers have also explored how system energetics impact flapping-wing performance. These studies provide foundational work for exploring the interactions between the drive and wing subsystems for real-world flapping-wing air vehicles. Researchers at the U.S. Army Research Laboratory and University of Maryland have conducted motor characterization together with load cell testing of their Robo Raven FWAV to characterize force production and motor performance under varying operational conditions [21,22]. Karpelson et al. used a strip theory aerodynamic model to explore the energetics of a very small hovering flapping-wing robotic insect and provide useful design insight including a rationale for battery sizing and flight endurance predictions [23]. Madangopal et al. studied the energetics of a larger flapping-wing air vehicle, with consideration to the mechanism's role in storing and releasing energy, and how this may impact system efficiency by reducing losses associated with wing deceleration at the end of each wing beat [24,25]. Doman et al. [26] studied the interactions of the propulsion system, drive mechanism, and flexible wing spars, revealing the system-level behavior of varying drive system parameters with a Lagrangian-approach to energy analysis. However, this study also highlights the challenges in validating model results with experimental testing, as the model becomes very complex despite the lack of an aerodynamic prediction component. All of these efforts reveal the importance of considering multiple vehicle subsystems together, due to interactive effects that dictate overall system performance.

In this work, a new model for the dynamics of a flapping wing is developed to improve motion predictions. Unlike previous work, drive motor constraints on torque and speed are included. This enables the wing motion to be computed by enforcing feasible motor bandwidth at each time step instead of simply prescribing plunge and twist motions. The deflection properties of the wing spars are captured in a dynamic model that is used to improve correlation between experimentally measured flapping dynamics and the modeled flapping dynamics. This approach addresses the need for a simultaneous consideration of motor bandwidth and flapping loads, since their interaction will have strong effects on overall flapping-wing performance. A modified strip theory modeling framework is used to predict aerodynamic loads and therefore is dependent on experimental validation data to tune the model parameters. To utilize accurate data during model calibration, an instrumentation suite that provides in-flight data is integrated into a flapping-wing aerial vehicle, known as Robo Raven II [21,27]. Using this data, model parameters are tuned to enforce feasible operational conditions. This approach addresses the problems associated with rigid mounting structures used in traditional laboratory experiments that restrict body responses to flapping. The model predictions show that flapping motions realized in flight testing significantly deviate from commanded sinusoidal motions due to interactions between the drive motors and wing loads. The new model for wing motions is coupled to the simplified aerodynamic model to determine the improvement in prediction accuracy that can be achieved by including subsystem interactions to enforce realistic wing motions with respect to drive motor constraints.

Aerodynamic Model for Flapping Wings Using Strip Theory

For the purposes of this study, a modeling approach with rapid solution throughput is desired to provide the opportunity to iteratively solve for feasible flapping motions at each time step. Therefore, the modified strip theory approach developed by DeLaurier is used as a model baseline [15]. Several corrections are applied to this model to improve the prediction quality in the present study which will be briefly summarized.

The model discretizes the wings into chordwise strips from the root to the tip, computes quasi-steady forces on each strip, and sums the forces to determine the total wing force. Wing chord is the product of the root chord defined as c0 and the chord shape c¯(ζ), defined in Eq. (1) as a function of normalized span location ζ. This wing design is derived from an experimental characterization and manufacturing sensitivity analysis that established a wing baseline with a favorable blend of lift and thrust production [28].
c¯(ζ)={1 for 0ζ<124ζ(1ζ)for 12ζ1
(1)
Stripwise normal force is computed as the sum of circulatory and added mass forces. Circulatory force is computed as shown in the below equation, where the flight speed and relative flow velocity and the quarter chord are U and V̂0.25c, semispan is y0, spanwise location is ζ, and the normal force coefficient is Cn
dNc=ρUV̂0.25c2Cn(y0ζ)cos(γ)cy0dζ
(2)
The normal force coefficient is computed using a modified three-dimensional Theodorsen function to account for finite-span unsteady vortex wake effects. The added mass force is computed as shown in the below equation
dNa=ρπc24(Uα˙14cθ¨)y0dζ
(3)
The chordwise forces are computed as the sum of cambered wing drag, leading-edge suction, and viscous friction drag as shown in the below equations
dDcamber=2πα0(α+θ¯)cos(γ)ρUV̂0.25c2cy0dζ
(4)
dTs=ηs2πα0(α+θ¯)cos(γ)ρUV̂0.25c2cy0dζ
(5)
dDf=CdfρVx22 cy0dζ
(6)
The baseline modified strip theory model used thus far determines if flow separation occurs following a helicopter analysis methodology that is unsuitable for the flapping-wing case now examined. Therefore, an enhanced dynamic stall criterion is adopted from the model by Kim et al. for large amplitude flapping [19] and based on experimental data collected by Scherer in characterizing oscillating airfoils [29]. The dynamic stall condition is defined in Eq. (7) with the dynamic stall correction factor defined in Eq. (8) as a function of both plunging and pitching
αdyn=ξmaxαstall
(7)
ξmax={1+|tan1[h˙cos(θθ¯a)/Vx]||αstall|+0.51θ|θ˙|c|θ˙|2Uαstall,αdyn2αstall2,αdyn<2αstall
(8)

The method of force application is dependent on flow condition at each time step. With attached flow, circulatory forces act at the quarter chord point, added mass acts at the half chord point, and chordwise forces act along the wing chord. In the case of dynamic stall, the leading-edge suction force shifts to the wing normal direction and the circulatory force moves to the 1/3 chord location. Finally, in separated flow conditions, all the chordwise forces vanish and all the normal forces act at the half chord point.

The differential forces in the normal and chordwise directions are resolved into lift and thrust forces at each time step based on the wing positions according to Eqs. (9) and (10), then the total wing forces are computed as the sum of each differential wing strip
dL=dNcosθ+dFxsinθ
(9)
dT=dFxcosθdNsinθ
(10)

The basic framework used for aerodynamic modeling and dynamic stall modeling has been presented here for completeness but the detailed description of the aerodynamic models chosen is beyond the scope of the present effort. However, highly detailed treatments of the modified strip theory aerodynamic model as well as the dynamic stall corrections applied are available from the respective authors, including several illustrations depicting the terms in the equations used and experimental validation studies [15,19].

Instrumentation Integration Into Flapping-Wing Flight Platform

Application of the presented strip theory approach requires flight testing to validate model parameters including vehicle orientation, airspeed, and wing tracking, corresponding to θ¯, U, and γ in aerodynamic model equations. To provide this information while eliminating the errors associated with rigid mounting of flapping-wing air vehicles, a free flight instrumentation suite is integrated into a flapping-wing air vehicle.

A high level diagram of the instrumentation suite is shown in Fig. 1. The system uses four parallel microcontrollers to handle computationally expensive tasks due to the need for processor availability to respond quickly to pilot inputs and the need for a stable timestamp across all the data channels collected by the system. The attitude and heading reference system (AHRS) uses sensor fusion algorithm to estimate vehicle attitude using a direction cosine matrix approach that eliminates sensor drift summarized in Ref. [30]. The analog to digital conversion (ADC) of current and voltage measurements is handled by a separate microcontroller to maximize sampling rates for improved power consumption estimates. This approach is necessary because the servo motors used have a control system with 300 Hz PWM current modulation instead of smooth analog adjustments. The main microcontroller synchronizes computed results from the AHRS and ADC microcontrollers with the remaining sensors, while also handling the radio receiver and servo drive for command and control of the vehicle. Finally, the results of all the synchronized sensor outputs and pilot commands are written to a MicroSD card by a dedicated processor due to the high data throughput required to prevent buffer overruns. As a standalone add-on, a GPS data logging module may be optionally included which consists of an antenna and data logging system. This module is removable to facilitate indoor flights in GPS-denied environments.

Fig. 1
Instrumentation suite high level functional diagram
Fig. 1
Instrumentation suite high level functional diagram
Close modal

Preliminary testing and calibration trials were conducted to ensure appropriate sensor measurement ranges and precision. All the selected sensor model numbers and calibration results are shown in Table 1. Due to the large number of required interconnections in the instrumentation suite, the sensors were integrated into a printed circuit board to minimize size, weight, and complexity. The final assembled printed circuit board weighs 7.5 g and is shown in Fig. 2.

Fig. 2
Custom PCB used in Robo Raven II flight tests
Fig. 2
Custom PCB used in Robo Raven II flight tests
Close modal
Table 1

Sensors used on Robo Raven II flapping-wing aerial vehicle

SensorModelVoltageRangePrecision
VoltageATMEGA328P5.00–5.5 V4.9 mV
CurrentACS7235.0±10 A20 mA
Optical encoderE2-5005.0N/A0.18 deg
AccelerometerADXL3453.3±16 gN/A
MagnetometerHMC5883L3.3±8 GN/A
GyroscopeITG-32003.3±2000 deg/sN/A
Diff. pressure (pitot)HSCMRRN001ND2A55.0±248.84 Pa0.2 m/s airspeed
GPSPA6H-MTK33393.3N/A3 m
Data loggingOpenLog3.3Up to 1 M Baud
SensorModelVoltageRangePrecision
VoltageATMEGA328P5.00–5.5 V4.9 mV
CurrentACS7235.0±10 A20 mA
Optical encoderE2-5005.0N/A0.18 deg
AccelerometerADXL3453.3±16 gN/A
MagnetometerHMC5883L3.3±8 GN/A
GyroscopeITG-32003.3±2000 deg/sN/A
Diff. pressure (pitot)HSCMRRN001ND2A55.0±248.84 Pa0.2 m/s airspeed
GPSPA6H-MTK33393.3N/A3 m
Data loggingOpenLog3.3Up to 1 M Baud

A dynamometer and power analyzer were used to conduct motor characterization following the approach outlined in prior work [21], resulting in models that map the voltage and current to motor operational characteristics. The torque model is shown in Eq. (11), and the power output model is shown in Eq. (12). The experimentally measured model parameters are listed in Table 2 for the selected drive servos

Table 2

Experimentally measured parameters for Futaba S9352HV servo

ParameterFutaba S9352HV
CA−4.688 W
CB4.646 W/A
CC0.677 W/V
CD−0.936 A
CE−6.877 V
CF0.919
CG−1.270 W/A2
CH−0.070 W/V2
ParameterFutaba S9352HV
CA−4.688 W
CB4.646 W/A
CC0.677 W/V
CD−0.936 A
CE−6.877 V
CF0.919
CG−1.270 W/A2
CH−0.070 W/V2
τ=KτI+Kτ0
(11)
Pout=CA+CBI+CCV+(I+CD)((V+CE)CF)+(I+CD)((I+CD)CG)+(V+CE)((V+CE)CH
(12)

A Robo Raven II flapping-wing aerial vehicle was equipped with the developed instrumentation suite and several test trials were conducted to establish the range of operational conditions in typical cruising flight. First, maneuvering dynamics were explored by executing a series of turns consisting of either a simple tail yaw motion, or a coordinated wing and tail maneuver that also included asymmetric flap amplitude changes to create a thrust differential. These tests were executed ten times each using a prerecorded maneuver sequence to ensure precise maneuver start and stop times at the same phase relative to the standard flapping motion, then the results were cycle-averaged to yield clear trends. Results of these tests are shown in Fig. 3.

Fig. 3
Maneuver test results
Fig. 3
Maneuver test results
Close modal

In addition to maneuvering flight, several trials of cruising flight were conducted. The experimental trials included four wings following the design template shown in Fig. 4 with dimensions listed in Table 3 and chord shape described by Eq. (1). A detailed description of this wing design including the deformation characteristics, performance properties, and a manufacturing sensitivity analysis is available in prior work [21,28]. For each wing design, several experimental trials were conducted with the objective of achieving a stable flight condition. This was achieved by making adjustments to the battery mounting location to shift the center of gravity of the vehicle. This changes the body orientation in flight and thus alters the flight speed and power consumption, which is expected to be a roughly U-shaped curve as is exhibited by many species of flying birds and bats [31].

Fig. 4
Generalized wing design template
Fig. 4
Generalized wing design template
Close modal
Table 3

Wings used in experimental validation

Wing designMass (g)Span (m)Max chord (m)Area (m2)
A12.70.5270.3110.113
B19.40.6220.3560.171
C21.00.6670.3940.210
D25.40.7620.4320.280
Wing designMass (g)Span (m)Max chord (m)Area (m2)
A12.70.5270.3110.113
B19.40.6220.3560.171
C21.00.6670.3940.210
D25.40.7620.4320.280

A snapshot of a typical flight test result is shown in Fig. 5, with 50 Hz sample rate and sample count plotted on the x-axis. The results show airspeed and altitude in the first row, with the flight starting at approximately sample number 500 where a spike in airspeed can be seen from the initial hand launch. The second row shows pilot control inputs during flight, the third row shows outputs from the AHRS orientation estimation, and the final row shows the voltage and current consumption. Summary results from all the testing trials are plotted in Figs. 57 and provided in Table 4. While the number of trials is still not sufficient to fully characterize the envelope of feasible performance, the trends emerging in the data set are extremely useful in detailing known feasible flight conditions. This insight is more useful than load cell testing that has been conducted in prior work [21,22] because it shows how flight dynamics impact vehicle functionality in real flight conditions.

Fig. 5
Cruising flight test results
Fig. 5
Cruising flight test results
Close modal
Fig. 6
Airspeed and inclination results from cruising flight tests
Fig. 6
Airspeed and inclination results from cruising flight tests
Close modal
Fig. 7
Power consumption in cruising flight at 4.0 Hz flapping
Fig. 7
Power consumption in cruising flight at 4.0 Hz flapping
Close modal
Table 4

Cruising flight summary results

ParameterWing AWing BWing CWing D
Climb rate (m/s)0.4890.7580.8510.395
Angle of attack (rad)0.8790.8390.3900.585
Airspeed (m/s)5.544.806.427.75
All-up weight (kg)0.3120.3290.3290.359
Flap amplitude (rad)1.1371.0621.0491.032
Mean current (A)2.542.852.933.26
ParameterWing AWing BWing CWing D
Climb rate (m/s)0.4890.7580.8510.395
Angle of attack (rad)0.8790.8390.3900.585
Airspeed (m/s)5.544.806.427.75
All-up weight (kg)0.3120.3290.3290.359
Flap amplitude (rad)1.1371.0621.0491.032
Mean current (A)2.542.852.933.26

Modeling of Wing Motion

A wing subsystem model is required that captures the important interactive effects between wing design, flapping motions, and force production. At present, models used to describe flapping wings such as [15] do not appropriately capture the interaction between flapping motions and wing compliance and loads. In particular, two key deficiencies exist. First, sinusoidal flapping is frequently assumed. This assumption does not consider the retarding effect of wing loads. When the wing loading is considered together with the motor capabilities, the realized motion is augmented based on the wing size and motor bandwidth and no longer may be assumed to be sinusoidal. The second deficiency in existing models is the static wing shape specification, which assumes a constant linear wing twist and does not accurately describe the behavior of a flexible wing that deforms in response to flapping and structural loads. Accounting for these deficiencies will improve the aerodynamic model quality by ensuring that real-world effects are appropriately considered in performance predictions.

In Fig. 8, the plunge motions are shown for each wing while flapping at 4.0 Hz. A counterintuitive effect takes place during pronation and supination where the larger wings are able to reach higher speeds and reduce the phase gap relative to the smaller wings. Following these portions of the flap cycle, the power stroke exhibits the expected behavior of larger wings falling behind the smaller wings due to larger loading arising from flapping a greater surface area. In addition, the motions measured for each wing are significantly different than the pure sinusoidal motion that is commanded. The dynamic deformation of the wing spar structure in response to loading will be explored in this section to provide a more accurate representation of flapping motions in the aerodynamic model.

Fig. 8
Wing angle tracking results
Fig. 8
Wing angle tracking results
Close modal

The wing deformation has two separate modes that depend on flapping velocity and wing stiffness. The first is the primary bending mode, which consists of the wing spars flexing perpendicular to the wing surface in response to loading. The second is a twisting mode relative to the primary leading-edge spar that occurs during stroke reversal. The bending mode controls drag which causes a torque load at the motor that scales linearly with flapping velocity and wing size and remains in phase with the flapping velocity. To illustrate this effect, the torque required by the motor to drive each wing design at varying steady-state angular velocity was recorded to establish the relationship between wing size and drag. A regression for each wing design is plotted in Fig. 9. The coefficients describing the relationship between wing size, flapping velocity, and torque required are listed in Table 5. Overlaid on the plot is the bandwidth for the chosen drive servos, which bounds a region of feasible operation for each wing size. Any operational conditions that fall below and to the left of this line are reachable by a given wing design.

Fig. 9
Torque required for steady plunge velocity
Fig. 9
Torque required for steady plunge velocity
Close modal
Table 5

Torque required for each wing design in steady-state plunge motion

WingArea (m2)Torque constant (N·m s/rad)
A0.16400.0595
B0.22130.1127
C0.26250.1439
D0.32900.2002
WingArea (m2)Torque constant (N·m s/rad)
A0.16400.0595
B0.22130.1127
C0.26250.1439
D0.32900.2002

In Fig. 10, the angular velocity is plotted for wing D at the 1.0 Hz and 4.0 Hz flapping rates. The angular velocities corresponding to commanded kinematics are plotted as solid lines, while the measured angular velocities are plotted as dots. The overlaid dashed lines are the torque limitations for the wing as shown in Fig. 9. Two important effects may be observed in this plot. First, the motor is able to exactly reach the commanded plunge profile as long as torque is maintained below the motor bandwidth. This is clear in the 1.0 Hz flapping condition where the commanded and actual motions are indistinguishable. However, as the commanded flapping motions become more demanding, the motor bandwidth constrains the reachable angular velocity.

Fig. 10
Comparison between actual and commanded angular velocity for wing D
Fig. 10
Comparison between actual and commanded angular velocity for wing D
Close modal

In Fig. 11, wing D angular velocity test results are shown for a range of flapping rates. At flapping rates below 1.5 Hz, there is a smooth sinusoidal velocity profile because the flapping motions are not sufficiently demanding to exceed motor bandwidth. Beyond 1.5 Hz, the motor bandwidth limit is reached, however, the angular velocity is not strictly limited by the theoretical maximum speed dictated by motor bandwidth. The angular velocity exhibits an overshoot effect that causes significantly higher flapping velocity for brief periods. High-speed videography was used to investigate the underlying physical reason for this effect. A series of snapshots from one flap cycle are shown in Fig. 12.

Fig. 11
Angular velocity for wing D across flap rates
Fig. 11
Angular velocity for wing D across flap rates
Close modal
Fig. 12
High-speed video captures wing D test
Fig. 12
High-speed video captures wing D test
Close modal

In the video snapshots, the flapping motion starts with the top row moving from left to right, then the bottom row moving from left to right. An upstroke is in progress in the top left frame. Moving along the top row, the upstroke is completed by the drive motor, and hence, the directly mounted primary wing spar reaches its apex in the fourth frame. On the fifth frame through the seventh frame, the wing rotates as stored elastic energy is released, which aids the beginning of the downstroke by reducing torque requirements temporarily. It is this rotational effect that augments the flapping velocity and causes the overshoot effect observed in higher flapping rates plotted in Fig. 11. Evidently, the twisting effect causes overshoot that depends on flapping frequency. This becomes important when the wing stroke reversal rate approaches the natural wing twisting dynamics. Once these two effects begin to overlap, the torque requirements are significantly reduced because the wing is exhibiting dynamic twisting and releasing stored elastic energy that aids the stroke reversal.

The twist dynamics are accounted for in the wing model by introducing several corrections to the nominal flapping profile. To setup these corrections, first a notional wing system model is constructed as shown in Fig. 13. In this model, the two axes are defined as the primary flapping axis f̂ and the wing twist axis t̂. The two generalized coordinates that describe wing motion are θf and θt, which are the angle of the primary spar at the leading edge of the wing relative to the horizontal and the wing twist angle relative to the plane that includes the leading-edge spar and the flapping axis. The energy system of the wing is modeled using torsion dampers attached to the flapping axis and the twist axis, a torsion spring attached to the twist axis, and a rotational inertia J with components in the flap and twist axes. Functionally, the effect of this model framework is that the wing plunging motion exhibits a dependence on the interaction between the flap rate, the wing size, and the twist stiffness.

Fig. 13
Two-axis wing flexibility model
Fig. 13
Two-axis wing flexibility model
Close modal
Using this model framework, the wing plunge rate is increased by the twist dynamics based on the interaction between wing natural frequency and flapping rate. As the flapping rate overlaps the wing natural frequency, the amount of twist will increase resulting in higher velocity peaks in plunge. To capture this effect, the wing natural frequency is described using an empirical relationship based on the wing design used. The primary structure resisting wing twist is the bending deformation of the spars in the chordwise direction; therefore, the twisting stiffness is modeled as proportional to cantilevered beam bending as described in Eq. (13). Since each wing uses constant spar sizes of the same carbon fiber material, the elastic modulus and second moment of area terms are lumped together with the constant Kt̂ to capture bending physics and the empirical constant of proportionality B together
Kt̂=3EIL3B
(13)
Combining the mass of each wing listed in Table 3 with the empirical stiffness relationship results in the wing natural frequency in the below equation
ωnt̂=Kt̂m
(14)
The twist is modeled as a damped vibratory system with damped natural frequency given by the below equation
ωdt̂=ωnt̂1ζ2
(15)
The solution for the homogeneous displacement response is shown in the below equation
θt̂(t)=θ0t̂eζωnt̂tcos(ωdt̂t)+θ˙0t̂+ζωnt̂θ0t̂ωdt̂eζωnt̂tsin(ωdt̂t)
(16)

The initial twist condition θ0t̂ is estimated using high-speed photography to determine the wing twist for each wing tested as a function of flapping rate. A snapshot of this testing is shown in Fig. 14. The damping ratio ζ is estimated in Eq. (17) as a function of flap rate per wing using the percent overshoot observed in Fig. 11 following the relationship in Eq. (17) [32], where the overshoot is measured by comparing the peak plunge rate achieved during pronation to the steady-state plunge rate limit after velocity ringing settles. The results of damping ratio calculations are compiled in Table 6 for each wing size tested

Fig. 14
High-speed photography used to characterize wing twist amplitude
Fig. 14
High-speed photography used to characterize wing twist amplitude
Close modal
Table 6

Wing twist damping ratio results

Wingθ˙max (rad/s)θ˙ss (rad/s)Percentage overshootζ
A9.499.470.200.96
B10.389.3710.800.58
C10.969.1220.180.45
D11.488.4136.500.31
Wingθ˙max (rad/s)θ˙ss (rad/s)Percentage overshootζ
A9.499.470.200.96
B10.389.3710.800.58
C10.969.1220.180.45
D11.488.4136.500.31
ζ=(lnPO100)2π2+(lnPO100)2
(17)

The predicted plunge velocity incorporating the wing twist augmentation model is plotted in Fig. 15 along with the commanded and experimentally measured angular velocity for wing D at 4.0 Hz flapping. The early inaccuracy in the model is due to a brief transient effect that settles quickly and results in a much closer fit to actual flapping motions by accounting for the interaction of the wing stiffness and twisting.

Fig. 15
Computed flap motion incorporating motor model
Fig. 15
Computed flap motion incorporating motor model
Close modal

Experimental and Modeling Results

The flapping-wing system is modeled by creating a linkage between the strip theory model and the motor torque model to perform feasibility checking at each time step as shown in Fig. 16. This is realized by computing the aerodynamic loads and flapping speeds associated with nominal flapping profile, mapping the loads to a torque and angular velocity bandwidth for the drive motors, and using feedback control to correct the flapping profile until feasible motor operation is achieved. If the loads are either too large to be driven at the current speed or less than what the motor is capable of, the plunge rate is reduced and forces recomputed until the solution converges to agreement between the two models, which prevents violation of the feasibility constraint. In this way, the modeling approach is a reflection of the digital control system used by the drive servo to minimize position error. This approach is a new strategy for modeling wing dynamics in flapping-wing air vehicles, since wing motions are predicted at each time step by considering the interactions between the motors and wings. In contrast, the traditional strategy has been to design a flapping mechanism which follows specified flapping kinematics and to represent the same kinematics exactly in the modeling without corrections due to motor loading conditions or wing dynamics.

Fig. 16
Strip theory model approach
Fig. 16
Strip theory model approach
Close modal

Each of the typical flight conditions corresponding to the four wing designs as listed in Table 4 was modeled using the aerodynamic code combined with the wing twist model and motor constraints to establish model capabilities in describing real flight conditions. The flight testing parameters were entered into the aerodynamic model to check the ability of the model to describe known flight-worthy conditions and to identify areas requiring tuning. A comparison of the flight test results to the modeled results both with and without the kinematic corrections is shown in Table 7.

Table 7

Comparison of model results to flight testing data


Flight testing

Coupled model (developed in this paper)

Kinematics only
ABCDABCDABCD
Avg. lift (N)2.503.083.193.312.514.034.756.42
Avg. thrust (N)1.161.151.372.191.162.883.957.15
Avg. power (W)13.615.115.917.612.918.320.022.113.029.839.970.9
Avg. torque (N·m)0.991.051.061.170.701.021.151.311.432.302.703.64
Max. plunge velocity (rad/s)14.6914.0413.6012.6116.6915.0212.9611.5620202020

Flight testing

Coupled model (developed in this paper)

Kinematics only
ABCDABCDABCD
Avg. lift (N)2.503.083.193.312.514.034.756.42
Avg. thrust (N)1.161.151.372.191.162.883.957.15
Avg. power (W)13.615.115.917.612.918.320.022.113.029.839.970.9
Avg. torque (N·m)0.991.051.061.170.701.021.151.311.432.302.703.64
Max. plunge velocity (rad/s)14.6914.0413.6012.6116.6915.0212.9611.5620202020

In the comparison, lift and thrust data are not populated for the flight test results since these are not directly measurable from the instrumentation suite. In the coupled model results, the inclusion of wing dynamics and a motor feasibility constraint causes significant reduction in the power predictions. This results in model conditions that much more closely track the behavior exhibited during the flight testing trials. For comparison, the model results that specify flapping kinematics without any modifications may be thought of as estimating the requirements to drive a given wing at exactly the nominal flapping kinematics. The corrected data present an estimate of the achieved performance given the interactions between the selected components. The major reason underlying the improvement in the couple model is the reduction in required motor bandwidth. The torque requirements necessary to achieve sinusoidal flapping kinematics are far beyond the capability of the motors selected, as shown by the discrepancy in the two models. In addition, the kinematics-only model prescribes the same sinusoidal flapping motion across each wing size, despite the significant differences in loading that must be overcome by the motors as wing area increases. With the modified wing motion model, the plunge velocity is reduced to a feasible condition that places the operation of the motor within feasible bounds and ensures that modeled kinematics are more representative of real flight conditions.

Conclusions

In this study, a coupled modeling approach has been developed that improves the prediction of wing motions during flapping by simultaneously considering drive motor constraints along with the aerodynamic loading conditions modeled using strip theory, rather than simply specifying trigonometric flapping kinematics. To provide data input for the aerodynamic model, an instrumentation suite has been developed and deployed on a flapping-wing air vehicle. Improved prediction of flapping motion is achieved because feasible operational conditions are enforced at each time step. The inclusion of these subsystem interactions reduces errors in the predicted wing motion by ensuring the motors operate in physically realistic conditions. Comparisons of the measured and predicted wing motion confirm the improvement relative to the idealized wing motions.

The wing deformation model treats the wing as a spring–mass–damper system with bending and torsional modes and uses high-speed photography to determine shape deformation in response to flapping. This method allows the aerodynamic model to include realistic flapping motions, however, it lacks a detailed linkage between unsteady aerodynamics and structural dynamics. Therefore, the proposed approach is suitable for generating performance predictions of the presented design family, but may suffer from inaccuracies when applied to new wing designs without conducting additional tuning of the wing dynamics model. More broadly, the approach presented here is useful in robotics applications that model the motions of dynamic end effectors. Often, motions are prescribed as a model initialization step, with no further regard to the interactive effects of time-varying forces. The approach offered here seeks to reduce errors in model outputs by considering how limited actuator capabilities together with dynamic interactions between forces and motions may alter the overall system behavior.

The approach presented in this paper provides improved wing motion predictions, resulting in reduced errors in the modified strip theory aerodynamic model. To extend the results to include measures of performance including endurance and range, a battery derating model is necessary to capture performance loss associated with battery discharge rate and state of charge. In addition, a limited set of test cases have been conducted so far to demonstrate some valid conditions for flight. To further improve the accuracy of the model, it will be necessary to generate a large matrix of test cases and statistically evaluate the performance variability. A final area for future work will be to leverage the present work into an interactive design tool with a dashboard-style layout similar to existing quadcopter and airplane tools.1

Acknowledgment

John Gerdes' participation on this research was sponsored by the U.S. Army Research Laboratory, Vehicle Technology Directorate. S. K. Gupta's and Hugh A. Bruck's participation on this research were supported by AFOSR Grant No. FA9550-12-1-0158, managed by Dr. Byung-Lip “Les” Lee.

References

1.
Pennycuick
,
C.
,
1990
, “
Predicting Wingbeat Frequency and Wavelength of Birds
,”
J. Exp. Biol.
,
150
(
1
), pp.
171
185
.
2.
Gerdes
,
J. W.
,
Wilkerson
,
S. A.
, and
Gupta
,
S. K.
,
2012
, “
A Review of Bird-Inspired Flapping Wing Miniature Air Vehicle Designs
,”
ASME J. Mech. Rob.
,
4
(
2
), p.
021003
.
3.
de Croon, G. C. H. E., Ruijsink, R., De Wagter, C., Perçin, M., and Remes, B. D. W., 2016,
The DelFly: Design, Aerodynamics, and Artificial Intelligence of a Flapping Wing Robot
, Springer, Dordrecht, The Netherlands.
4.
Zdunich
,
P.
,
Bilyk
,
D.
,
MacMaster
,
M.
,
Loewen
,
D.
,
DeLaurier
,
J.
,
Kombluh
,
R.
,
Low
,
T.
,
Stanford
,
S.
, and
Holeman
,
D.
,
2007
, “
Development and Testing of the Mentor Flapping-Wing Micro Air Vehicle
,”
J. Aircr.
,
44
(
5
), pp.
1701
1711
.
5.
Keennon
,
M.
,
Klingebiel
,
K.
,
Won
,
H.
, and
Andriukov
,
A.
,
2012
, “
Development of the Nano Hummingbird: A Tailless Flapping Wing Micro Air Vehicle
,”
AIAA
Paper No. 2012-0588.
6.
Pornsin-Sirirak
,
T.
,
Tai
,
Y.
,
Ho
,
C.
, and
Keennon
,
M.
,
2001
, “
Microbat: A Palm-Sized Electrically Powered Ornithopter
,”
NASA/JPL
Workshop on Biomorphic Robotics
, Pasadena, CA, Aug. 14–17.
7.
Send
,
W.
,
Fischer
,
M.
,
Jebens
,
K.
,
Mugrauer
,
R.
,
Nagarathinam
,
A.
, and
Scharstein
,
F.
,
2012
, “
Artificial Hinged-Wing Bird With Active Torsion and Partially Linear Kinematics
,” 28th Congress of the International Council of the Aeronautical Sciences (
ICAS
), Brisbane, Australia, Sept. 23–28.
8.
Mueller
,
D.
,
Bruck
,
H. A.
, and
Gupta
,
S. K.
,
2009
, “
Measurement of Thrust and Lift Forces Associated With Drag of Compliant Flapping Wing for Micro Air Vehicles Using a New Test Stand Design
,”
Exp. Mech.
,
50
(6), pp.
725
735
.
9.
Hubel
,
T. Y.
, and
Tropea
,
C.
,
2009
, “
Experimental Investigation of a Flapping Wing Model
,”
Exp. Fluids
,
46
(
5
), pp.
945
961
.
10.
Shkarayev
,
S.
, and
Maniar
,
G.
,
2013
, “
Experimental and Computational Modeling of the Kinematics and Aerodynamics of Flapping Wing
,”
J. Aircr.
,
50
(
6
), pp.
1734
1747
.
11.
Shkarayev
,
S.
, and
Silin
,
D.
,
2010
, “
Aerodynamics of Cambered Membrane Flapping Wings
,”
AIAA
Paper No. 2010-58.
12.
Song
,
A.
,
Tian
,
X.
,
Israeli
,
E.
,
Galvao
,
R.
,
Bishop
,
K.
,
Swartz
,
S.
, and
Breuer
,
K.
,
2008
, “
Aeromechanics of Membrane Wings With Implications for Animal Flight
,”
AIAA J.
,
46
(
8
), pp.
2096
2106
.
13.
Caetano
,
J.
,
Percin
,
M.
,
Oudheusden
,
B. V.
,
Remes
,
B.
,
Wagter
,
C. D.
,
Croon
,
G. D.
, and
Visser
,
C. D.
,
2015
, “
Error Analysis and Assessment of Unsteady Forces Acting on a Flapping Wing Micro Air Vehicle: Free Flight Versus Wind-Tunnel Experimental Methods
,”
Bioinspiration Biomimetics
,
10
(
5
), p.
056004
.
14.
Grauer
,
J. A.
, and
Hubbard
,
J. E.
,
2008
, “
Development of a Sensor Suite for a Flapping-Wing UAV Platform
,”
AIAA
Paper No. 2008-224.
15.
DeLaurier
,
J.
,
1993
, “
An Aerodynamic Model for Flapping-Wing Flight
,”
Aeronaut. J.
,
97
(
964
), pp.
125
130
.
16.
Jones
,
R. T.
,
1940
, “
The Unsteady Lift of a Wing of Finite Aspect Ratio
,” NACA Technical Report, Report No. NACA-TR-681.
17.
Mazaheri
,
K.
, and
Ebrahimi
,
A.
,
2012
, “
Performance Analysis of a Flapping-Wing Vehicle Based on Experimental Aerodynamic Data
,”
J. Aerosp. Eng.
,
25
(
1
), pp.
45
50
.
18.
Hunsaker
,
D. F.
, and
Phillips
,
W. F.
,
2015
, “
Propulsion Theory of Flapping Airfoils, Comparison With Computational Fluid Dynamics
,”
AIAA
Paper No. 2015-0257.
19.
Kim
,
D.-K.
,
Lee
,
J.-S.
, and
Han
,
J.-H.
,
2011
, “
Improved Aerodynamic Model for Efficient Analysis of Flapping-Wing Flight
,”
AIAA J.
,
49
(
4
), pp.
868
872
.
20.
Smith
,
M. J. C.
,
Wilkin
,
P. J.
, and
Williams
,
M. H.
,
1996
, “
The Advantages of an Unsteady Panel Method in Modelling the Aerodynamic Forces on Rigid Flapping Wings
,”
J. Exp. Biol.
,
199
(
5
), pp.
1073
1083
.
21.
Gerdes
,
J. W.
,
Bruck
,
H. A.
, and
Gupta
,
S. K.
,
2015
, “
A Systematic Exploration of Wing Size on Flapping Wing Air Vehicle Performance
,”
ASME
Paper No. DETC2015-47316.
22.
Gerdes
,
J. W.
,
Roberts
,
L.
,
Barnett
,
E.
,
Kempny
,
J.
,
Perez-Rosado
,
A.
,
Bruck
,
H. A.
, and
Gupta
,
S. K.
,
2013
, “
Wing Performance Characterization for Flapping Wing Air Vehicles
,”
ASME
Paper No. DETC2013-12479.
23.
Karpelson
,
M.
,
Whitney
,
J.
,
Wei
,
G.-Y.
, and
Wood
,
R.
,
2010
, “
Energetics of Flapping-Wing Robotic Insects: Towards Autonomous Hovering Flight
,”
2010 IEEE/RSJ International Conference on Intelligent Robots and Systems
(
IROS
), Taipei, Taiwan, Oct. 18–22, pp.
1630
1637
.
24.
Madangopal
,
R.
,
Khan
,
Z.
, and
Agrawal
,
S.
,
2004
, “
Energetics Based Design of Small Flapping Wing Air Vehicles
,”
IEEE International Conference on Robotics and Automation
, (
ICRA
), New Orleans, LA, Apr. 26–May 1, pp.
2367
2372
.
25.
Madangopal
,
R.
,
Khan
,
Z.
, and
Agrawal
,
S.
,
2005
, “
Biologically Inspired Design of Small Flapping Wing Air Vehicles Using Four-Bar Mechanisms and Quasi-Steady Aerodynamics
,”
ASME J. Mech. Des.
,
127
(
4
), pp.
809
817
.
26.
Doman
,
D. B.
,
Tang
,
C. P.
, and
Regisford
,
S.
, “
Modeling Interactions Between Flexible Flapping Wing Spars, Mechanisms, and Drive Motors
,”
J. Guid. Control Dyn.
,
34
(
5
), pp.
1457
1473
.
27.
Gerdes
,
J. W.
,
Holness
,
A.
,
Perez-Rosado
,
A.
,
Roberts
,
L.
,
Greisinger
,
A.
,
Barnett
,
E.
,
Kempny
,
J.
,
Lingam
,
D.
,
Yeh
,
C.-H.
,
Bruck
,
H. A.
, and
Gupta
,
S. K.
,
2014
, “
Robo Raven: A Flapping-Wing Air Vehicle With Highly Compliant and Independently Controlled Wings
,”
Soft Rob.
,
1
(
4
), pp.
275
288
.
28.
Gerdes
,
J. W.
,
Cellon
,
K.
,
Bruck
,
H. A.
, and
Gupta
,
S. K.
,
2013
, “
Characterization of the Mechanics of Compliant Wing Designs for Flapping-Wing Miniature Air Vehicles
,”
Exp. Mech.
,
53
(
9
), pp.
1561
1571
.
29.
Scherer
,
J. O.
,
1968
, “
Experimental and Theoretical Investigation of Large Amplitude Oscillation Foil Propulsion Systems
,” Hydronautics, Inc., Laurel, MD, Report No.
TR-662-1-F
.
30.
Bartz
,
P.
,
2016
, “
Building an AHRS Using the SparkFun “9DOF Razor IMU” or “9DOF Sensor Stick”
,” GitHub, Inc., San Francisco, CA, accessed Aug. 10, 2016, https://github.com/ptrbrtz/razor-9dof-ahrs/wiki/tutorial#mathematical-background-and-firmware-internals
31.
Rayner
,
J. M. V.
,
2001
, “
Mathematical Modelling of the Avian Flight Power Curve
,”
Math. Methods Appl. Sci.
,
24
(
17–18
), pp.
1485
1514
.
32.
Ogata
,
K.
,
2009
,
Modern Control Engineering
, 3rd ed.,
Pearson
,
Upper Saddle River, NJ
.

Nomenclature

     
  • B =

    empirical stiffness constant

  •  
  • c =

    chord

  •  
  • C =

    motor output power empirical constants

  •  
  • c0 =

    chord at wing root

  •  
  • Cn =

    coefficient of normal force

  •  
  • Cdf =

    coefficient of friction drag

  •  
  • c¯ =

    chord shape function

  •  
  • D =

    drag

  •  
  • E =

    Young's modulus

  •  
  • Fx =

    chordwise force

  •  
  • h =

    wing plunge angle along flapping axis

  •  
  • I =

    current

  •  
  • Kt̂ =

    torsional stiffness constant

  •  
  • Kτ =

    motor torque empirical constants

  •  
  • L =

    lift

  •  
  • m =

    wing mass

  •  
  • N =

    normal force

  •  
  • Pout =

    motor output power

  •  
  • PO =

    percent overshoot

  •  
  • t =

    time

  •  
  • T =

    thrust

  •  
  • U =

    airspeed

  •  
  • V =

    voltage

  •  
  • Vx =

    horizontal velocity on wing section

  •  
  • V̂0.25c =

    resultant velocity at quarter-chord

  •  
  • y0 =

    semispan

  •  
  • α =

    relative angle of attack at 3/4-chord

  •  
  • αdyn =

    dynamic stall angle

  •  
  • αstall =

    static stall angle

  •  
  • α0 =

    zero-lift angle

  •  
  • α =

    angle of attack at 3/4-chord due to unsteady effects

  •  
  • γ =

    dihedral

  •  
  • ζ =

    nondimensional span coordinate

  •  
  • ηs =

    suction efficiency

  •  
  • θ =

    chord pitch relative to U

  •  
  • θt̂ =

    wing twist deflection

  •  
  • θ0t̂ =

    wing twist initial condition

  •  
  • θ¯ =

    mean pitch angle relative to U

  •  
  • θ¯a =

    inclination of flapping axis relative to U

  •  
  • ξmax =

    stall angle coefficient

  •  
  • ρ =

    air density

  •  
  • τ =

    torque

  •  
  • ωdt̂ =

    torsional damped natural frequency

  •  
  • ωnt̂ =

    torsional natural frequency