IEEE/ASME TRANSACTIONS ON MECHATRONICS, VOL. 18, NO.

2, APRIL 2013 787

A Simulation and Design Tool for a Passive Rotation

Flapping Wing Mechanism
Veaceslav Arabagi, Lindsey Hines, and Metin Sitti, Senior Member, IEEE

Abstract—This paper develops a numerical simulation tool for wing flapping only, while allowing for passive rotation of the
designing passive rotation flapping wing mechanisms. The simu- wing.
lation tool includes a quasi-static model of the piezoelectric ac- Wood et al. proposed a flapping wing platform capable
tuator bending, transmission kinematics, and the small Reynolds
number aerodynamic forces governing wing dynamics. To validate of fully controlled wing trajectory, and thus body forces and
the developed tool, two single-wing systems with distinct resonant torques [9]. However, that robot design relies on a difficult to
frequencies are manufactured and characterized. Comparison to manufacture differential mechanism, driving both the wing’s
experimental results reveals that, although discrepancies exist, the leading and trailing edges, and a complex high-bandwidth
simulation is able to predict general trends of wing kinematics and controller [10]. However, recent theoretical and experimental
lift behavior as functions of frequency, thus, being useful as a design
tool. Finally, the complex models, ranging from actuator deflection results of Bergou et al. suggest that at least some of the wing
to wing aerodynamics presented in this paper, allow analysis of rotation in insects is passively assisted by inertial and aero-
the complete system revealing insight into several wing trajectory dynamic forces [1]. Furthermore, the flapping wing regime is
control methodologies, and potentially serving as a design and op- known to possess several self-stabilizing damping properties for
timization tool for future flapping wing robots. translation and rotation maneuvers due to the reciprocal nature
Index Terms—Flapping wing, passive wing pitch reversal, sim- of the wing beat [11], [12], and upright maintenance with the aid
ulation and design tool, wing lift force modulation. of sail type devices [13]; thus, suggesting the possibility of an
underactuated wing drive mechanism. Based on this idea, a pas-
I. INTRODUCTION sive wing rotation approach has been undertaken by Wood, with
a single actuator driving the flapping motion of the wings while
S in natural flyers, miniature flapping wing vehicles are
A expected to present a significant advancement in agility
from their fixed and rotary wing counterparts. Inspiration and
allowing for passive rotation, with the amplitude of rotation
limited by mechanical stoppers [14]. This mechanism results
in large impact forces on the wing imparted at the end of each
basic concepts of flapping flight are taken from studying the
stroke, which could decrease its lifetime and excite undesirable
kinematics of dragonflies [1], butterflies [2], hawkmoths [3],
structural vibration modes.
bats [4], and flies [5]. The low Reynolds number operation
This paper is based upon an aerial platform design featur-
regime of small-scale fliers is dominated by nonlinear aerody-
ing completely passive wing pitch reversal. Without mechanical
namics and temporal wing–wake interactions [6], [7] that are
stops and suspended from an axis close to its leading edge by a
the underlying cause for their, unaccounted by conventional
spring/damper system, the motion of the wing is governed only
aerodynamics, large lift force production capacity [8], respon-
by the dynamics of the wing, aerodynamic forces, and the elas-
sible for the acrobatic maneuverability of these insects. One
ticity of the system. Furthermore, this design not only simplifies
major milestone to overcome in creating a robotic flapping
mechanical complexity, but also presents opportunity for the
wing platform capable of controlled hover is the generation
eventual introduction of wing angle of attack control via a vari-
of proper and controlled wing motion. Over time, two differ-
able spring stiffness at the wing rotation joint, potentially using
ent approaches to this problem emerged: 1) active control of
smart, morphing materials. A dynamic simulation that includes
both the flapping and rotation angles, and 2) active control of
a compilation of small Reynolds number aerodynamic forces,
coupled with dynamics of the driving piezoelectric bending ac-
tuator and kinematics of the four-bar transmission, serves as the
basis for analysis of the system’s performance. This model is
Manuscript received February 8, 2011; revised September 2, 2011 and
October 27, 2011; accepted December 21, 2011. Date of publication February employed to simulate the behavior of two single-flapping-wing
20, 2012; date of current version January 10, 2013. Recommended by Technical systems, one resonating at 30 Hz and the other at 75 Hz, and the
Editor S. Fatikow. The work of L. Hines was supported by the National Defense results are compared to experimental measurements. Although
Science and Engineering Graduate Fellowship.
V. Arabagi was with the Department of Mechanical Engineering, Carnegie the theoretical model is not able to predict the exact lift force
Mellon University, Pittsburgh, PA 15213 USA. He is now with the Department magnitudes, wing rotation dynamics are predicted accurately
of Cardiovascular Surgery, Children’s Hospital Boston, Boston, MA 02115 USA with experimental flapping angles as input. In addition, the sim-
(e-mail: veaceslav.arabagi@childrens.harvard.edu).
M. Sitti is with the Department of Mechanical Engineering, Carnegie Mellon ulation is able to capture behavioral trends very well, proving
University, Pittsburgh, PA 15213 USA (e-mail: msitti@andrew.cmu.edu). useful as a simulation and design tool. Finally, three methods
L. Hines is with the Robotics Institute, Carnegie Mellon University, to control the lift force of a passively rotating flapping wing
Pittsburgh, PA 15213 USA (e-mail: lhines@cmu.edu).
Color versions of one or more of the figures in this paper are available online mechanism are proposed: modulating the flapping frequency,
at http://ieeexplore.ieee.org. the amplitude of the actuator oscillation, and the torsional stiff-
Digital Object Identifier 10.1109/TMECH.2012.2185707 ness of the wing rotation joint.
1083-4435/$31.00 © 2012 IEEE
788 IEEE/ASME TRANSACTIONS ON MECHATRONICS, VOL. 18, NO. 2, APRIL 2013

Fig. 2. Displacement (left axis) and torque transmission curves (right axis)
for 75-Hz resonance single-flapping-wing system. τ o u t is the torque output of
Fig. 1. Schematic of the complete single-wing system, including the actuator,
the four-bar. The lengths of these four-bar’s links are {L0, L1, L2, L3} = {13.5,
slider crank, four-bar, and wing.
12.1, 10, 5} mm.

II. SYSTEM DYNAMICS MODELING

bending, a slider crank mechanism is implemented to transfer
A numerical simulation of the actuator-transmission-wing the actuator output forces into the transmission. The four-bar
dynamics is developed in MATLAB to serve as a design and joints are modeled as ideal pin joints, with the exception of the
optimization tool. The simulation includes three main compo- joint at the output link of the four-bar, designated by θ in Fig. 1,
nents: the actuator, a transmission mechanism, and the wing. which has a stiffness value assigned to it, as it undergoes the
largest range of motion and hence can significantly affect system
A. Actuator dynamics. Based on the four-bar transmission kinematics mod-
Piezoelectric actuators are selected to power the robotic flap- eled in [23], the output angle and torque transmission function
per due to their high efficiencies, power densities of up to of the four-bar is
400 W/kg [15], extensive theoretical modeling [16], [17], easy
θ = h(linksi , conf ig, δ) (2)
incorporation into the composite manufacturing process, and
their use with some success by other groups, such as the Berke-  drive = f (linksi , conf ig, Fin )
M (3)
ley MFI [18], [19] and the Harvard MicroFly [20]. Although
piezoelectric materials also have several limitations, such as a where M  drive is the moment that is applied to the wing about
need for a high-voltage power supply, small strain (less than its flapping axis (E3 axis from Fig. 3), linksi are the lengths
0.1% for PZT-5H), and high brittleness, work has been done of the 4 four-bar links, conf ig is the configuration of the four-
to make their application feasible in miniature aerial vehicles bar, i.e., its position through the stroke, and Fin is the force
as presented in [21]. Reciprocating wing drive is achieved with input from the actuator, as portrayed in Fig. 1. These general
a bending unimorph piezoelectric actuator, whose schematic is equations are sufficient to describe any transmission system
portrayed in Fig. 1. The modeled actuator consists of an ac- transforming linear input displacement into angular output mo-
tive bending portion and a rigid extension that increases its tion. In the case of a piezoelectric actuator, Fin = kact (Δ − δ)
stroke. Classical lamination theory was used to model the ac- and represents how the actuator stiffness enters into the theo-
tuator quasi-statically, alike presented in [17]. The structural retical model. Numerical optimization is employed to achieve a
elasticity of the unimorph actuator is introduced into the model particular transmission ratio and symmetry over the actuator dc-
as a spring of stiffness kact attached to the end of a rigid ac- bias voltage, yielding transmission curves portrayed in Fig. 2.
tuator that bends under an applied electric field, as shown in Generally, higher torque transmission ratios imply a higher res-
Fig. 1. The inertial properties of the actuator have been simpli- onance frequency if everything else is kept constant. Intuitively,
fied to the dynamics of a point mass positioned on the tip of this arises from the fact that higher transmission ratios impart
the actuator. The effective mass of the actuator, meff , calculated smaller wing inertial forces on the actuator and hence the latter
based on an energy approach presented in [22], was determined senses a smaller effective wing inertia.
to be 15 mg for the 30-Hz resonant system. Its kinetic energy,
used for the Lagrangian expression, is Ta = 12 meff δ̇ 2 , where δ C. Wing Aerodynamics
is the displacement of the slider crank end. The free, no load,
displacement of the actuator tip is modeled as Passive rotation of the wing is achieved by attaching the wing
to a driver spar with an elastic joint, thus, allowing the elas-
Δ = g(E, actGeom) (1) tic/damping properties of the flexure material to influence wing
motion. In fact, the interaction of these elastic and damping
where E is the electric field in the PZT layer and actGeom is
torques with the wing’s inertia and aerodynamic forces is essen-
the geometry of the actuator.
tial to the generation of proper flapping wing trajectory. Fig. 3
portrays the geometric setup of the simulation model and the
B. Transmission
sets of Euler angles for coordinate transformations.
A four-bar-type transmission is employed to amplify the small Throughout its motion, the wing is affected by aerodynamic
tip displacements of the piezoelectric actuator. Due to actuator forces stemming from various small Reynolds number effects.
ARABAGI et al.: SIMULATION AND DESIGN TOOL FOR A PASSIVE ROTATION FLAPPING WING MECHANISM 789

ρπ
dFair = − (θ̈cosϕ − θ̇ϕ̇sinϕ)rc(r)2 dr
4
 
ρπ c(r)
− ϕ̈ zRA − c(r)2 dr (7)
4 2
where ρ is the air density, r is the radial position of a wing strip
from the wing’s flapping axis, c(r) is the chord length of the
particular strip, θ and ϕ are the flapping and rotation angles,
respectively (portrayed in Fig. 3), α is the angle of attack of
the wing, or (π/2 − ϕ), and zRA is the location of the rotation
axis measured from the wing’s leading edge. The translational
force component is obtained by vector addition of mutually
Fig. 3. Schematics of the passive flapping wing setup. Coordinate sets repre- orthogonal lift and drag forces. U is taken to be the velocity in
sent transformations established by θ (flapping angle) and ϕ (rotation angle).
The coordinate systems are shifted for clarity, while in simulation they are all the E2 direction of each wing strip’s midchord point
centered at the point labeled “coordinate sets origin.” The E coordinate system
  = rθ̇ + c(r)
is attached to the wing. U (r) = vm c · E 2 ϕ̇cosϕ (8)
2
where vm c is the midchord velocity of each wing division. The
use of midchord velocity in the translational lift force expres-
sion is justified by the fact that the incoming air stream velocity
is linearly distributed along each chord of the wing; hence, the
best approximation to a constant velocity field, which was used
to obtain the translational lift coefficients in [8], is the mean ve-
locity along each wing chord, observed at the midchord point.
Cl , Cd , and Crot are the translational lift, drag, and rotational
force coefficients, respectively, where the translational coeffi-
cient expressions were tabulated experimentally by Sane and
Dickinson for a Reynolds number (Re) of 136 [25]

Fig. 4. (a) Wing SolidWorks model and (b) wing profile illustration of aero- Cl (α) = 0.225 + 1.58 sin(2.13α − 7.2) (9)
dynamic force due to rotation. The wing center of mass and rotation axis are
illustrated. Cd (α) = 1.92 − 1.55 cos(2.04α − 9.82) (10)
where α is in degrees. The position of the rotation axis has been
Overall, the total force on the wing is
set at ≈1/4 chord length from the leading edge of the wing, at
Ftot = Ft + Frot + Fair + Fwc (4) its aerodynamic center (according to thin-plate airfoil theory).
The Reynolds number for our system is Re = 2800, calculated
where Ft is the translational force, Frot is the force due to wing by averaging the Re numbers obtained for each wing strip, is
rotation, Fair is the force due to added air mass, and Fwc is the similar, in aerodynamic sense, to RoboFly’s one, Re = 136,
force created by wake capture. Translational force estimates are justifying the use of Sane and Dickinson’s findings as means of
quasi-steady approximations adapted from thin airfoil theory. approximating the aerodynamic forces on the wing. The value
Frot is the force generated by air resistance to wing rotation as for the rotational drag coefficient Crot was taken to be 2, as
portrayed in Fig. 4(b). The added air mass forces are imparted this is the theoretical result of rotational drag on a flat plate
on the wing by the surrounding air and are generated due to subjected to normal flow [26], [27]. To account for the radial
unsteady wing motion [24]. This is also the only force that aids variability of wing velocity, the aforementioned equations need
the rotation of the wing during stroke reversal, since the inertia to be integrated numerically to obtain the aerodynamic forces
of the air trapped in the vicinity of the wing tends to continue on the entire wing.
its linear motion, thereby imparting a torsional moment on the
wing [25]. The wake capture effect is a nonsteady phenomenon
D. Equations of Motion of the Passive Rotation Wing System
occurring when the wing traverses vortices and air circulation
generated by its motion prior to the direction change. Since a The wing has a carbon fiber leading edge spar with an ex-
closed form expression for this force cannot be determined, it is truding vein for chordwise stiffness and a Kapton membrane
not included in the dynamical model. Thus, the expressions for as the main airfoil (see Fig. 4). The inertial properties of the
aerodynamic forces acting on chordwise wing strips, dr, take wing required for the dynamic equations are obtained directly
the form from the SolidWorks model. In addition, the three carbon fiber
1
dFt = ρU 2 c(r)[Cl2 (α) + Cd2 (α)]1/2 dr (5) links of the four-bar and the slider crank mechanism were mod-
2 eled as inertial elements. The dynamic equations are formulated

1 in Lagrangian form, with generalized coordinates of θ and ϕ
dFrot = − Crot ρ|ϕ̇|ϕ̇ z  |z  |dz  dr (6)
2 c(r ) being the flapping and rotation angles, respectively. Thus, we
790 IEEE/ASME TRANSACTIONS ON MECHATRONICS, VOL. 18, NO. 2, APRIL 2013

formulate the Lagrangian as follows: moment due to the aerodynamic forces, explicitly defined as
1 1 
L=T −V = mv · v + J ω·ω   aero =
M i (α)
Fi × β (16)
2 2
i
1 1
+ mL i vL i · vL i + JL i ωL i · ωL i
2 2 where Fi ’s are the translational and rotational aerodynamic
forces and βi ’s are the respective positions of their centers of
1 1
+ meff δ̇ 2 − krot ϕ2 (11) pressure from the wing rotation axis, explicitly defined as
2 2
 
where m is the wing mass augmented with the added air mass 0.82
[as described in (19)], v is the velocity at the wing’s center |βt (α)| = |α| + 0.05 c(r) (17)
π
of gravity (CG), J is the wing’s inertia matrix taken about the  
wing’s center of mass, krot is the torsional spring stiffness at the span z dFrot
|βrot | =  (18)
rotation axis, mL i , JL i , vL i , and ωL i are the masses, rotational span dFrot
inertias, linear and angular velocities of the four-bar links 1–3,
respectively, meff is the effective mass of the actuator, δ̇ is the where βt defines the position of the translational lift center of
horizontal velocity of the tip of the actuator (scalar), and ω  is pressure as a function of angle of attack for each wing chord
the angular velocity of the wing defined as strip c(r) and is obtained from experimental results of Dickson
et al. [28]. βrot is the effective moment arm of the rotational
ω   + ϕ̇E
 = θ̇E   = ϕ̇E
  + θ̇ sin(ϕ)E
  + θ̇ cos(ϕ)E
  (12)
3 1 1 2 3 force distribution of Fig. 4(b).
expressed in terms of basis vectors of the coordinate frame at- Given that there is no analytical solution for the added air
tached to the wing. Furthermore, the expressions for angular mass force on a wing planform moving in 3-D fashion, dFair
and linear velocity of the four-bar links are complex in closed of (7) does not accurately describe the moments exerted by
form and, hence, have been substituted by numerically com- the added air mass on the flapping and rotation angles. As an
puted trajectories parametrized by the output flapping angle approximation, the effect of added air mass forces was imple-
θ (i.e., v2 (t) = f (θ(t)) and hence dvdt
2 (t)
= θ̇(t) ∂ ∂θ (t)
f
). Defin- mented in the form of virtual mass that augmented the physical
 ˙ mass of the wing, similar to the inertial implementation in [26]
ing the conventions f ≡ ∂f
and f ≡ , and substituting
df
∂ θ (t) dt and [29]. Stemming from the concept that the added mass effect
all knowns into the Lagrangian L, we obtain the following can be estimated as dm = π4ρ c(r)2 dr [30], the expression for
equation: m emerges as follows:
L = 1/2(−krot ϕ2 + sin(ϕ)2 Jy y θ̇2 
πρ
m = mw + c(r)2 dr (19)
+ ϕ̇(cos(ϕ)Jxz θ̇+Jxx ϕ̇)+cos(ϕ)θ̇(cos(ϕ)Jz z θ̇ + Jxz ϕ̇) span 4
2
+ 1/2m((2RCG + β 2 cos(2ϕ))θ̇2 + 4RCG βcos(ϕ)θ̇ϕ̇ where mw is the physical mass of the wing. This approximation

3 
2 allows a simplified treatment of the added air mass effects that

+ 2β 2 ϕ̇2 ) + mL i θ̇2 xL i 2 + JL i θ̇2 θL2i is sufficient for the development of a simulation tool able to
i=1 i=1 predict general dynamics and lift force trends. Although this
2 2 
2 mass augmentation approach is used in the differential equations
+ JL 3 θ̇ + meff θ̇ δ ) (13) defining the system dynamics, the expression for added air mass
where RCG and βCG are the radial and vertical positions, re- force of (7) is used to estimate the aerodynamic lift generated
spectively, of the center of mass from the origin of the coordinate by this effect in all lift plots presented in this paper.
systems, as portrayed in Fig. 4, Jii are the components of the The previously presented Lagrange equations (14) and (15)
wing’s inertia matrix, xL i is the CG position of link i of the are fully general, governing the dynamics of a passive rotation
four-bar, and θL i is the orientation of the links (note that the flapping wing driven by a torque M  drive around the flapping
four-bar motion is planar and equations are simplified to reflect axis. In our case of the driving mechanism consisting of a piezo-
2-D motion). Note that for a thin wing, the inertial components electric actuator coupled to a four-bar transmission, M  drive is
Jy z and Jxy are extremely small and have been omitted from given by (3). Thus, substituting all the known quantities into (14)
the equations. Thus, the Lagrange equations describing the wing and (15) and simplifying, we obtain the nonlinear differential
flapping and rotation dynamics are as follows: equations governing the flapping and rotation angles that can be
  found in the Appendix, due to their length and complexity.
d ∂L ∂L   − dϕ̇
 aero · E
− =M 1 (14)
dt ∂ ϕ̇ ∂ϕ
  III. MANUFACTURING AND EXPERIMENTAL
d ∂L ∂L  drive + M  aero ) · E
 3 − Dθ̇ (15) PARAMETER DETERMINATION
− = (M
dt ∂ θ̇ ∂θ In order to best quantify the theoretical simulation’s perfor-
where d is the rotational damping coefficient, D is the flap- mance, two experimental systems consisting of an actuator, a
ping damping coefficient arising from the transmission flexure four-bar and wing were manufactured, having different reso-
damping, M drive is the driving flapping torque, and M
 aero is the nant frequencies of 30 and 75 Hz. Both systems feature the
ARABAGI et al.: SIMULATION AND DESIGN TOOL FOR A PASSIVE ROTATION FLAPPING WING MECHANISM 791

were recorded with a laser scan micrometer (Keyence LS-3100,

Woodcliff Lake, NJ), and, being second-order responses, al-
lowed the determination of stiffness and damping parameters
from the dominating frequency and envelope of the signal. These
test were performed under vacuum such as to eliminate aerody-
namic damping, which would constitute the main component of
damping otherwise.

D. Transmission Compliance Characterization

Given that the slider crank and four-bar flexures are sources
Fig. 5. Overall dimensions of the manufactured wing and actuator for the
of unwanted compliance, the stiffness and damping properties
(a) 30-Hz resonant system and (b) 75-Hz system. The actuator mount is 3-D of the transmission are measured experimentally. Compliance
printed from VisiJet HR 200 material, and assembled with Cyanoacrylate glue. of the slider crank is measured via an indentation test at the tip
of the first transmission link, with the actuator rigidly fixed. The
same construction components, differing only in their dimen- recorded tip stiffness is transformed via the transmission ratio
sions and lengths. to stiffness of the slider crank and modeled in simulation in
series with the actuator tip stiffness kact . Throughout our exper-
A. Component Manufacturing iments, slider crank stiffness ranges as low as ksc = 300 N/m
and as high as ksc = 1650 N/m. The large span of stiffness val-
All the parts comprising the flapping wing system were man- ues is attributed to the easy misalignment of the slider crank
ufactured using the “Smart Composite Microstructures” (SCM) links, resulting in flexure buckling. For all simulations in this
methodology and manually assembled afterward [9]. This last paper, the 30-Hz resonance flapping system was simulated with
manual step generally introduces some misalignment and uncer- a nominal slider crank stiffness of ksc = 1400 N/m, while the
tainty in the final assembly, which tends to adversely affect the 75-Hz system had a value of ksc = 300 N/m, per experimental
dynamics of the system, generally in the form of unwanted com- measurements.
pliance. The transmission members are formed of two preim- The damping coefficient of the four-bar was measured by im-
pregnated, unidirectional M60J ultrahigh modulus carbon fiber parting an impulse on wing-driving mechanism, with the wing
layers of 60-μm thickness (Toray) sandwiching a 12.7-μm layer rotated 90◦ , i.e., α = 0◦ , so as to eliminate aerodynamic drag
of Kapton (DuPont Kapton 50HN) and cured together. The wing, from the measurements. The simulated impulse response was
whose Kapton layer serves as a membrane, is manufactured from then matched to experimentally recorded vibrations via altering
thinner Kapton of merely 6.3 μm (DuPont Kapton 25HN), as the flapping damping parameter D from (15), converging to a
the minimization of the wing weight is of major importance. final value of D = 100 μN·mm·s, employed in simulation.
The actuator’s active layer consists of two layers of carbon
fiber bonded to a 125-μm-thick layer of PZT-5H (Piezo Sys- IV. EXPERIMENTAL SETUP AND WING KINEMATICS
tems), while the passive extension has an additional layer of
s-glass, increasing stiffness. Some dimensions of the actuator Experimental lift was measured by mounting the single-wing
and the wing used in the experimental flapping wing systems assembly to a load cell through a mechanical amplification link-
are portrayed in Fig. 5. age to increase the sensed lift forces, and the low-friction, single-
degree-of-freedom (DOF) pivot isolates the vertical forces gen-
B. Actuator Experimental Characterization erated by the wing [see Fig. 6(b)]. A 25-g load cell (Transducer
Techniques) for the 75-Hz resonant prototype (30-g load cell
Since the piezoelectric actuator is of key importance in defin- for the 30-Hz system) was attached to the main shaft via a
ing system dynamics, experiments were performed to obtain the contact connection, such that no torques in the shaft would be
actual actuator blocking force, tip stiffness, and material prop- transmitted into the load cell. The actuator, wing, and four-bar
erties. Using laminar plate theory, the actuator bending stiffness mechanism were mounted directly to a 3-D printed (Invision
with a force applied at its tip, labeled kact in Fig. 1, replicating HR 3D Printer) plastic mount that was bolted onto an aluminum
the operational scenario, was calculated at 160 N/m. Through shaft, portrayed in Fig. 6(a). Given the large inertia of the bal-
indentation experiments, the actual actuator stiffness was found ance system, its resonant frequency was determined at 13 Hz.
to be 169 N/m, measured by the ratio of force over given tip This is sufficiently distant from the operating frequencies of
displacements. 30 and 75 Hz to ensure lift measurements were unexaggerated;
however, functioning as a low-pass filter, measurements are re-
C. Wing Rotational Flexure Characterization stricted to mean lift force only.
The dynamic parameters of the rotational flexure were exper- During each test, the input voltage to the actuator, a sinusoidal
imentally measured after manufacturing, to account for any im- wave specified by the peak-to-peak voltage Vpp and a constant
perfections. The stiffness and damping properties of each flexure dc offset Vdc . Wing kinematics were captured using a strobe light
were determined by applying an impulse impact to it, after an and a high-definition camcorder (Sony HDR-SR11). Video was
additional known mass was added. The resulting displacements taken from the top view of the wing, with the camera placed
792 IEEE/ASME TRANSACTIONS ON MECHATRONICS, VOL. 18, NO. 2, APRIL 2013

Fig. 6. Photographs of the experimental lift acquisition setup with main com-
ponents labeled. (a) Actuator-four-bar-wing experimental prototype mounted
onto a 3-D printed plastic part. (b) Lift balance mechanical amplification
device.
Fig. 7. (a) Overlayed experimental and simulated wing kinematics, and
(b) trajectories of simulated lift forces of the 30-Hz resonance single-wing
in an offset plane parallel to the stroke plane. Three points on system. System parameters: V p p = 180 V, V d c = 30 V, frequency = 30 Hz,
the wing were tracked, based on edge contrast, from which the k sc = 1400 N/m, k ro t = 20 mN·mm, and d = 10 μN·mm·s.
wing trajectory was calculated.
with the input of a sinusoidal signal driving the actuator, enables
A. Wing Kinematics Comparison analysis of the wing driving mechanism. This characterization
Given the previously described setup to obtain kinematics, of the simulation performance is done for both the 30- and
a quick comparison is performed to ensure that the numerical 75-Hz resonant systems, in order to give an idea about the
simulation behaves similar to the system it models: experimen- model’s shortcomings in simulating systems of different reso-
tal kinematic wing data are compared to predicted results in nances.
Fig. 7(a) for the 30-Hz resonant system. The kinematics feature
a discrepancy of ≈10◦ peak to peak in flapping angles, propa- A. Analysis of the 30-Hz Resonant System
gating in a discrepancy of ≈20◦ peak to peak in rotation, likely
Since a frequency response characterizes the behavior of a
due to the asymmetry of the torque transmission curve of the
dynamic system well, this methodology is applied to analyze
physical four-bar transmission. The simulated lift force compo-
the performance of the numerical model. For the 30-Hz reso-
nents shown in Fig. 7(b) portray that translational lift constitutes
nant experimental system, the frequency response was obtained
most of the total lift force with a small contribution from added
by recording the wing kinematics and load cell output while
air mass lift. The rotational lift component is very low, instead
sweeping the frequency of the input voltage signal at discrete
the effect of this force is manifested mostly through the torque
steps of 3 Hz. This relatively large value of the frequency step-
it imparts on the wing, thereby having a great impact on shaping
size is acceptable due to the low-quality factor (Q ≈ 2) of a
the kinematics of wing rotation and influencing total lift force
flapping wing system, ensuring smoothness of the frequency re-
in that manner.
sponse curve and, thus, no important missed dynamics. In order
to validate the quasi-steady aerodynamic models and underlying
V. RESULTS AND DISCUSSION
rotation dynamics, the leading edge of the simulation wing is
Characterization of the developed theoretical model of a driven based on experimental kinematics. The simulation takes
single-wing-flapping system is attempted through comparison as input for the flapping angle a sinusoid signal with the ampli-
to experimentally obtained wing kinematics and mean lift force tude obtained from measured wing kinematics. The kinematics
measurements. The full dynamic model can be analyzed by sep- plot portrays the two fundamental resonant frequencies of the
arating it into two subsystems: 1) the actuator/transmission drive system manifested by peaks in the flapping amplitude, occur-
and wing-flapping dynamics; 2) the aerodynamic force models ring around 30 and 36 Hz, as shown in Fig. 8. The two resonant
and wing rotation dynamics. The performance of the second points of the system correspond to an optimal balance between
subsystem can be quantified by using the experimentally mea- the internal stresses of the actuator, inertial forces, and nonlin-
sured wing flapping trajectory as a displacement source in the ear aerodynamic damping effects acting on the wing. Simulated
simulation of wing rotation only. Running the full simulation, wing rotation closely resembles the experimental observations
ARABAGI et al.: SIMULATION AND DESIGN TOOL FOR A PASSIVE ROTATION FLAPPING WING MECHANISM 793

Fig. 9. Frequency sweep obtained with the complete system simulation, in-
cluding the actuator drive and wing dynamics. System parameters: V p p = 180 V,
V d c = 30 V, k sc = 1400 N/m, k ro t = 20 mN·mm, and d = 10 μN·mm·s. Curves
labeled with a slider crank stiffness are obtained with simulated systems having
critical values k sc = 300 and 1650 N/m.

Fig. 8. Frequency sweep of the 30-Hz resonance system simulating only wing
rotation, with flapping angle set by a sinusoid wave with experimentally obtained
amplitudes. (a) Experimental and simulated wing kinematics, and (b) lift forces.
System parameters: V p p = 180 V, V d c = 30 V, k ro t = 20 mN·mm, and d =
10 μN·mm·s.

except for high frequencies; this is thought to occur due to the

excitement of second resonance mode and out-of-phase rota-
tion, addressed in detail in Section V-C. The lift force plots of
Fig. 8(b) portray very similar magnitudes as the experimental
mean lift curve. Translational lift is the dominant aerodynamic
force throughout this frequency range, with the added air mass
lift comprising half of this amount. Similar to the kinematics, the
simulated lift force curve diverges from experiments at high fre-
quencies, suggesting the dominance of unmodeled aerodynamic
effects present in that operating regime.
Simulating the complete flapping wing system along with
actuator and transmission when driven by a sinusoid signal,
we obtain the plots presented in Fig. 9. Both the simulated
kinematics and total lift force are comparable in magnitude to
their experimental measurements for frequencies below the sec-
ond resonance peak. The model predicts well the occurrence
of the first resonant peak while overpredicting the occurrence
of the second peak at 50 Hz, although being highly damped.
Fig. 10. Frequency sweep of the 75-Hz resonant system simulating only wing
The simulated frequency sweep is performed for a larger range rotation, with flapping angle set by a sinusoid wave with an experimentally
of frequency values, while experimental measurements were obtained amplitude. (a) Experimental and simulated wing kinematics, and
stopped earlier due to extreme wing rotations after the second (b) lift forces. System parameters: V p p = 150 V, V d c = 50 V, k ro t =
37.5 mN·mm, and d = 13 μN·mm·s.
resonance, destructive to the prototype. Furthermore, the curves
labeled with a slider crank stiffness value, ksc , represent total
lift force curves obtained from simulation runs of single-wing B. Analysis of the 75-Hz Resonant System
systems featuring bounding values for the corresponding values To obtain the 75-Hz resonant system frequency response, lift
of slider crank stiffness, namely ksc = 300 and 1650 N/m. The was measured with a balance system (see Fig. 6) when the ac-
change in the magnitude of the total lift curve as well as the tuator was driven via a continuous frequency sweep. Fig. 10(a)
observed shift in system resonance illustrates the system’s sen- portrays good matching between the experimental wing rota-
sitivity to sources of compliance in the wing drive mechanism. tion angles and the simulation modeling only wing rotation.
794 IEEE/ASME TRANSACTIONS ON MECHATRONICS, VOL. 18, NO. 2, APRIL 2013

Fig. 12. Frequency spectrum of the system impulse response. The motion of
the first four-bar link was recorded after an impulse has been applied. The first
and second system resonant peaks are indicated by arrows. Note that the 50-
and 100-Hz peaks of the experimental system are due to electrical noise in the
environment and do not reflect physical resonant frequencies.

system featuring a stiff slider crank. The other limiting value of

stiffness ksc = 300 N/m is the nominal magnitude for this exper-
Fig. 11. Frequency sweep obtained with the complete simulation. System
parameters: V p p = 150 V, V d c = 50 V, k sc = 300 N/m, k ro t = 37.5 mN·mm, imental prototype. Again, a change in this parameter propagates
and d = 13 μN·mm·s. The lift bounding curve corresponding to a slider crank into a resonant frequency increase by more than 10 Hz, having
stiffness of k sc = 300 N/m coincides with the nominal frequency response, great impact on the behavior of the system.
while the other bounding curve corresponds to k sc = 1650 N/m.

C. Discussion
The experimental flapping kinematics portray two resonance The presented comparison between experimental measure-
peaks spaced close to each other, similar to the low resonant ments and both wing rotation only and the complete system
frequency system. The simulated lift forces in Fig. 10(b) are simulation effectively illustrates the capabilities and limitations
slight underestimates from experimental results. The largest of the theoretical model of the single-flapping-wing system.
lift contribution comes from added air mass lift, with transla- The overall good agreement to experimental results of both
tional lift being half of its value. This result arises from the fact kinematics and lift forces produced by the simulation of wing
that the high-frequency low-amplitude operation regime of the rotation only suggests that wing inertial properties, flexure stiff-
75-Hz system maximizes the accelerations of the wing, and thus ness, and damping properties, as well as the rough, quasi-steady
the air in its vicinity, i.e., added air mass, while large flapping aerodynamic models are largely correct. The complete simu-
amplitudes and translational wing motion dominate the behavior lation, encompassing the actuator drive, yields less acceptable
of the 30-Hz resonant system. The large increase in simulated results, overpredicting lift force, especially for the 75-Hz reso-
lift at 95 Hz seen in Fig. 10(b), not present in the experimental nant frequency system. Given that most of this discrepancy is
response, arises from the increase in flapping amplitude at that attributed to the underprediction of flapping amplitude and that
frequency. This occurs due to the excitement of the second res- wing rotation dynamics are correct, as established earlier, the
onance mode in the experimental system, predicted poorly by underlying cause of the problem is traced to the flapping mo-
the simulated system, due to unmodeled aerodynamic effects, tion of the wing. The decrease in wing flapping displacement
discussed in detail in Section V-C. Furthermore, after examin- could arise from actuator depolarization, inconsistencies in the
ing the phase shift between the flapping and rotation sinusoids, manual alignment and assembly of the four-bar transmission, as
a shift of 81◦ is observed in simulation and one of 138◦ in well as from spanwise aerodynamic moments not encompassed
experiment, suggesting that only the experimental system has by quasi-static 2-D approximations.
begun the transition to out-of-phase rotation, thus, reducing the One of the limitations of the developed theoretical simulation
produced lift. is its inability to model accurately the behavior of the system for
The full simulation of the 75-Hz resonance system, portrayed frequencies higher than the first resonance. These discrepancies
in Fig. 11, overall portrays much worse agreement of simula- are thought to be products of nonlinear aerodynamic effects
tion to experiment. The twofold discrepancy in total lift force caused by large stroke and rotation amplitudes. As a test to
produced is attributed primarily to the greatly underpredicted this hypothesis, the 75-Hz resonant system’s impulse response
flapping amplitude in simulation, ≈20◦ peak to peak. The larger is recorded, featuring small flapping and rotation stroke am-
amplitudes also cause an eightfold increase in translational lift plitudes (<5◦ ) and thus rendering the aforementioned dynamic
from the results of Fig. 10(b), rendering it the dominant lift effects not important. Comparing the frequency components
component, alike in the 30-Hz resonance system. Similar to of the experimental and simulated impulse responses obtained
previous results, the simulation predicts well the first resonance with the fast Fourier transform algorithm, shown in Fig. 12, we
mode at 75 Hz, while overestimating the second one at 115 Hz, observe overlapping peaks at 78 and 129 Hz, suggesting that
while experimentally it is observed at 95 Hz. The curve of the inertial, force transmission, and stiffness properties of sys-
ksc = 1650 N/m, of Fig. 11, represents simulated total lift of a tem components were modeled adequately in the numerical
ARABAGI et al.: SIMULATION AND DESIGN TOOL FOR A PASSIVE ROTATION FLAPPING WING MECHANISM 795

simulation. Furthermore, the second resonance peak is de- imum lift is exaggerated, its occurrence at 10–20 Hz higher
creased from 129 to 95 Hz in actuated flapping experiments, than the first resonance peak is well predicted by the simulation
suggesting a large effect of added air mass and potentially other and enforced by experimental measurements. Hence, although
unmodeled forces on wing dynamics. The two resonance peaks an optimization of system parameters based on the numerical
should not be viewed as the flapping and rotation resonances model might yield overestimated lift forces, the physical opti-
of the wing, but rather as different modes of vibration of the mal point will lie in the close vicinity of the one predicted, as it
system. Given that these resonances are reflected as peaks in is defined by trends and general behavior, captured sufficiently
flapping amplitude, they portray an optimally balanced state by the numerical model. Furthermore, given that the theoreti-
between internal actuator stresses and aerodynamic forces im- cal model captures well the aerodynamic effects present in the
parted on the wing. The first resonance achieves this balance large amplitude flapping regime, it should be used as an opti-
via approximately symmetric wing rotation, with wing pitch mization routine for flapping wing systems designed to operate
reversal occurring roughly at the time of stroke reversal (the in those conditions. Moreover, given that large flapping ampli-
exact timing of rotation is a complex function of many differ- tudes are featured by most natural flapping wing insects, this
ent parameters), thus resulting in a considerable amount of lift feature should be sought in future system designs, rendering the
force. The second resonant peak maximizes flapping amplitude simulation tool extremely useful. The optimization of such sys-
because the drag on the wing significantly decreases, thus, al- tems should be done in accordance with both resonance modes
lowing larger deflections of the actuator. The reduction of drag of the system (such as the ones portrayed in Fig. 12), with their
occurs because the wing begins to rotate out of phase with flap- proximity dictated by optimality of lift force production at first
ping, such that the flapping and rotation curves of Fig. 7 are 180◦ resonance.
offset in phase. This onset of second resonance, when the center
of mass of the wing performs very little translational motion, is
the cause for the sharp decrease in experimental lift in the two VI. PROPOSED WING CONTROL METHODOLOGY
systems shortly after the first resonant frequency. The reason The passive wing rotation design allows proper wing motion
for the simulation’s overprediction of the second resonant fre- without the weight of additional actuators but, incorporated into
quency is attributed to unmodeled transient aerodynamic forces, a vehicle, results in an overall underactuated system. Deliberate
such as wake capture [6], [25]. These dynamic effects generated underactuated design is not uncommon (see for example un-
by the wing’s past trajectory create impulse-like loading on it deractuated mechanical grippers [31]), but can complicate the
during each stroke, thus, exciting the second resonance mode control of the final platform. A variety of control techniques
much earlier than predicted. exist for certain classes of nonlinear underactuated systems, in-
It is generally desirable to engineer the system such as to space cluding combinations of hybrid and switching control, passivity
the two resonance modes far apart in the frequency domain. based methods, coordinate transformations, and basic lineariza-
One may attempt to do so by utilizing a stiff wing rotational tion [32]. Adaptive control and model learning techniques have
joint, such as not to excite the second mode; however, this is been employed in restricted DOF scenarios [33]; however, as of
suboptimal from the performance perspective, since the wing yet, no overarching control theory for these systems has been
will rotate insufficiently in its first resonance mode, which is the developed, making the success of any particular system unsure.
main operating point. Furthermore, stiffening either the wing Nevertheless, control over the produced wing lift force is a vital
rotation flexure or the actuator has the effect of shifting both requirement for the incorporation of a flapping wing mechanism
peaks of Fig. 12 to higher frequencies, as well as increasing into any robotic platform, since future vehicle design could rely
their separation. The proximity of the two resonance peaks in primarily on this method for the production of body torques. We
an optimized lift system is inherent in the passive wing rotation propose three methods of modulating wing lift: controlling the
system at hand, where the inertia of the wing is significant and flapping frequency, the waveform amplitude driving the actua-
plays a major role in pitch reversal. By reducing the overall tor, and the torsional stiffness of the passive rotation joint. The
inertia/mass of the wing and employing a correspondingly stiff effectiveness of these methods is investigated experimentally,
wing rotation joint, it is possible to increase the frequency of resulting in data for Figs. 13 and 14. Note that these experi-
the second resonance much higher, thus, reducing interference ments were performed with a different single-wing prototype
and broadening the frequency range of maximum lift. than the rest of the experiments in this paper, having a resonant
Finally, although the model generally overpredicts total lift frequency of ≈70 Hz; however, the mechanical construction
and overestimates the system’s second resonance frequency, was identical to the prototypes discussed earlier.
its usefulness lies in capturing trends in wing kinematics as a Figs. 9 and 11 portray the system response to changing flap-
function of frequency, amplitude, spring stiffness, etc. Thus, ping frequency. From the control standpoint, overall lift in-
the model can be used as a design tool to obtain dimensions creases linearly until the resonant frequency, after which point
and configuration of the actuator, four-bar transmission, and the behavior diverges based on system construction. Further-
wing that would ensure wing lift is maximized at a particular more, given that system efficiency is maximized at resonance,
frequency. The first resonance peak is predicted accurately in only slight deviations from the main operating frequency are
simulation, but, since it is sensitive to slider crank stiffness, envisioned, rendering frequency a poor control input.
it should be designed in accordance with the bounding values Alternatively, varying the amplitude of the driving waveform
of that parameter. Although the predicted magnitude of max- of the actuator modulates the flapping stroke of the wing, thus,
796 IEEE/ASME TRANSACTIONS ON MECHATRONICS, VOL. 18, NO. 2, APRIL 2013

Fig. 15. Experimental lift surface. Represents variation of produced lift with
Fig. 13. Simulated and experimental system response to a sweep in the input input voltage and torsional spring stiffness. Experimental parameters: V d c =
waveform amplitude. Simulation parameters: V d c = 50 V, frequency = 70 Hz, 50 V, frequency = 70 Hz, and d = 5 μN·mm·s.
k ro t = 50 mN·mm, and d = 5 μN·mm·s. Experimental lift data are averaged
over three datapoints.
ture of the passive flapping wing allows slow actuation of the
rotational joint stiffness in flight. Five flexure joints of vary-
ing thickness were manufactured and tested in the experimental
system in order to quantify the effect of flexure stiffness. Exper-
imental and simulated results are portrayed in Fig. 14 and illus-
trate the expected correlation of decreasing lift with increasing
flexure joint stiffness. The increase in flapping amplitude with
increasing joint stiffness is attributed to the shifting of both res-
onance modes (as illustrated in Fig. 12 to higher frequencies and
an increase in Q, a fact experimentally confirmed. Control over
the stiffness of the passive rotation joint is envisioned via an
actively stiffness changing or smart material, such as dielectric
elastomers [35] or ionic polymer metal composite (IPMC) actu-
ators [36] actuators, with design, integration being future works.
Generally, it should be noted that, although each presented
control method has limited application and effect, a combination
of these mechanisms could be employed in a dual wing robot.
Indeed, control over rotational flexure stiffness and actuator in-
Fig. 14. Simulated system response to a sweep in the rotational spring stiff- put voltage allows production of any lift force on the surface
ness. Simulation parameters: V p p = 150 V, V d c = 50 V and frequency = 70 Hz, of Fig. 15, thus, maximizing the range of achievable lift forces.
and d = 5 μN·mm·s. Experimental lift data are averaged over three datapoints.
The lift force surface resembles a plane, suggesting a simple
linear controller implementation. Asymmetrical lift and, there-
affecting the lift force. The simulated system response to chang- fore, body torques could be achieved by altering the rotational
ing peak-to-peak driving amplitude, overlayed with experimen- stiffness on only one passive rotation wing joint of a dual wing
tal data, is portrayed in Fig. 13. A good correlation between aerial platform via a morphing polymer. Furthermore, due to a
experimental and simulated data is observed, although simu- reduced number of actuation inputs into the system, some of the
lated lift is consistently higher than in experiment. The linear six DOF will likely be coupled and uncontrolled, requiring either
behavior of the overall lift force and large span of produced additional actuation mechanisms or combining control inputs to
forces renders this control scenario very promising. Intuitively, achieve full position/orientation of the robot in free space. It
an increase in the lift force could always be achieved by in- should also be noted that, although being very advantageous in
creasing the flapping amplitude; yet, this approach is limited by some respects, the passive pitch reversal wing, proposed in this
electrical depoling and stress induced crack propagation of the paper, would render the aerial platform less agile and mobile
piezoelectric actuator. than its active control counterparts. However, the mechanical
The third envisioned method of controlling lift is modulating simplicity, reduced complexity of the controller, and naturally
the stiffness of the rotational flexure that suspends the wing. stable behavior of wing rotation, that benefit a passive rotation
Unlike the high-bandwidth control architecture required for the design, could potentially bring us closer to the goal of liftoff
Micromechanical Flying Insect [34], the inherently stable na- and stable hovering.
ARABAGI et al.: SIMULATION AND DESIGN TOOL FOR A PASSIVE ROTATION FLAPPING WING MECHANISM 797

VII. CONCLUSION involve inertial effects of the driving subsystem

It is of great importance for any design project to have a simu- θ̈ (mw RCG βCG cos(ϕ) + Jxz cos(ϕ))
lation tool to predict the performance of the final system. Given
that our goal is the development of a miniature, controlled, flap- − θ˙2 cos(ϕ)sin(ϕ)(mβCG
2
+ Jy y − Jz z )
ping flight robot, a full simulation of its subsystem—a single- 2   − dϕ̇E
 aero · E  
+ ϕ̈(mβCG + Jxx ) + ϕkrot = M 1 1 (20)
wing design—could greatly speed up the design process, as
well as enable optimization of its components. The numeri- where all the constant definitions are presented in the main body
cal simulation presented in this paper includes a model of the of this paper.
bending piezoelectric actuator, transmission kinematics, and the The equation governing the motion of θ is more complex as it
passively rotating wing aerodynamics. Two experimental pro- involves numerous inertial terms and parameterizations in terms
totypes featuring different resonance frequencies were manu- of θ. We need to remind ourselves that according to the notation
factured from carbon fiber composite material, enabling the discussed in the text, we have the following:
comparison of simulation to experimental results. Overall, there
was broad agreement between the flapping and rotation trajec- ∂f ∂2 f df d2 f
f ≡ , f  ≡ , f˙ ≡ , f¨ ≡ 2 . (21)
tory trends between the numerical simulation and experiment. ∂θ(t) ∂θ(t)2 dt dt
The major discrepancy between the measured and predicted lift
forces is attributed to manufacturing/material imperfections of Then, the equation of motion for θ becomes
the actuator and/or four-bar transmission, and unmodeled aero-   
β2
dynamic forces. However, the model is able to predict the main θ̈ sin2 (ϕ)Jy y + cos2 (ϕ)Jz z + RCG2
+ CG m
2
resonance peak accurately along with general trends of wing
kinematics and lift behavior as functions of frequency. Thus, 1 3
  
given its limitations in predicting the magnitude of produced − β 2 cos(2ϕ)m + mL i (xL2i + yL2i ) + meff δ 2
2 i=1
lift force, the simulation can be used as a design tool for opti-
mization of wing shape, actuator and four-bar geometry, as well 
2 

as a future dual-wing flapping robotic platform. Although the + JL i θL2i + JL 3 + ϕ̈cos(ϕ)(mRCG βCG + Jxz )
presented experimental prototype features a lift-to-weight ratio i=1
of ≈1/6 (measured on a robot employing the same construc- − ϕ̇2 sin(ϕ)(mRCG βCG + Jxz ) + θ̇ϕ̇sin(2ϕ)(mβCG
2
tion), scaled down variants of this design have demonstrated
liftoff [14]. Thus, given the direct scalability of the SCM man- 
3

ufacturing technique and similarity of Reynolds number, the + Jy y − Jz z ) + mL i θ̇2 (xL i xL i + yL i yL i )
i=1
presented here theoretical model and analysis is valid and use-
ful for scaled down robotic platforms capable of liftoff. 
2

Compared to other wing rotation ideologies, a completely + JL i θL i θL i θ̇2 + meff δ  δ  θ̇2
passive rotation design benefits from mechanical simplicity, i=1

passive stability of wing trajectory, and low controller band-  drive + M

= (M  aero ) · E
 3 − Dϕ̇. (22)
width. However, the high dependence of wing dynamics on
aerodynamic forces renders the system as lacking in robustness The coupled equations for θ(t) and ϕ(t) are solved numeri-
and consistent performance, requiring an adequate controller to cally with MATLAB’s ode15s solver. The obtained results are
compensate for these design issues. Control over the amplitude presented in the main body of the paper.
of the actuator driving waveform, as well as, the potential use of
smart materials in the wing rotation joint, is envisioned to enable ACKNOWLEDGMENT
indirect control over the trajectory of the wing. The proposed
The authors would like to thank R. Smith for his work on the
passive design is by no means seen as a replacement for active
mean lift measuring setup, as well as prototype manufacturing
wing control mechanisms, whose ability to generate arbitrary
and testing.
forces and torques are essential for execution of their tasks, but
merely as an alternative flapping platform that may be realizable
on a shorter time scale for specific applications. REFERENCES
[1] A. Bergou, S. Xu, and Z. Wang, “Passive wing pitch reversal in insect
flight,” J. Fluid Mech., vol. 591, pp. 321–337, 2007.
APPENDIX [2] R. B. Srygley, A. L. R. Thomas, “Unconventional lift-generating mecha-
nism in free-flying butterflies,” Nature, vol. 420, pp. 660–664, 2002.
FULL EQUATIONS OF MOTION FOR THE [3] T. L. Hedrick and T. L. Daniel, “Flight control in the hawkmoth manduca
PASSIVE ROTATION WING sexta: The inverse problem of hovering,” J. Exp. Biol., vol. 209, pp. 3114–
3130, 2006.
Given the definition for the Lagrangian in (13) and perform- [4] F. T. Muijres, L. C. Johansson, R. Barfield, M. Wolf, G. R. Spedding, and
A. Hedenstrom, “Leading-edge vortex improves lift in slow-flying bats,”
ing the differentiations in (14) and (15), the full equations of Science, vol. 319, pp. 1250–1253, 2008.
motion of wing motion are obtained. The equation governing [5] S. Fry, R. Sayaman, and M. H. Dickinson, “The aerodynamics of hovering
the rotation angle ϕ is the simplest of the two as it does not flight in drosophila,” J. Exp. Biol., vol. 208, pp. 2303–2318, 2005.
798 IEEE/ASME TRANSACTIONS ON MECHATRONICS, VOL. 18, NO. 2, APRIL 2013

[6] J. Birch and M. H. Dickinson, “The influence of wing-wake interactions [33] N. Perez-Arancibia, J. Whitney, and R. Wood, “Lift force control of
on the production of aerodynamic forces in flapping flight,” J. Exp. Biol., flapping-wing microrobots using adaptive feedforward cancelation sche-
vol. 206, pp. 2267–2272, 2003. mes,” IEEE/ASME Trans. Mechatronics, [Online]. Available:
[7] M. Hamamoto, Y. Ohta, K. Hara, and T. Hisada, “A fundamental study of http://ieeexplore.ieee.org, DOI: 10.1109/TMECH.2011.2163317.
wing actuation for a 6-in-wingspan flapping microaerial vehicle,” IEEE [34] X. Deng, L. Schenato, and S. Sastry, “Flapping flight for biomimetic
Trans. Robot., vol. 26, no. 2, pp. 244–255, Apr. 2010. robotic insects: Part II: Flight control design,” IEEE Trans. Robot., vol. 22,
[8] S. P. Sane and M. H. Dickinson, “The aerodynamic effects of wing rotation no. 4, pp. 789–803, Aug. 2006.
and a revised quasi-steady model of flapping flight,” J. Exp. Biol., vol. 205, [35] P. Lotz, M. Matysek, and H. Schlaak, “Fabrication and application of
pp. 1087–1096, 2002. miniaturized dielectric elastomer stack actuators,” IEEE/ASME Trans.
[9] R. J. Wood, S. Avadhanula, M. Menon, and R. S. Fearing, “Microrobotics Mechatronics, vol. 16, no. 1, pp. 58–66, Feb. 2011.
using composite materials: The micromechanical flying insect thorax,” in [36] C. Zheng and T. Xiaobo, “A control-oriented and physics-based model for
Proc. IEEE Int. Conf. Robot. Autom., 2003, pp. 1842–1849. ionic polymer–metal composite actuators,” IEEE/ASME Trans. Mecha-
[10] S. Avadhanula and R. S. Fearing, “Flexure design rules for carbon fiber tronics, vol. 13, no. 5, pp. 519–529, Oct. 2008.
microrobotic mechanisms,” in Proc. Int. Conf. Robot. Autom., 2005,
pp. 1579–1584.
[11] B. Cheng and X. Deng, “Translational and rotational damping of flapping
flight and its dynamics and stability at hovering,” IEEE Trans. Robot., Veaceslav Arabagi received the B.Sc. degree from
vol. 27, no. 5, pp. 1–16, Oct. 2011. the University of California, Berkeley, in 2006,
[12] T. L. Hedrick, B. Cheng, and X. Deng, “Wingbeat time and the scaling of and the M.Sc. and Ph.D. degrees from Carnegie
passive rotational damping in flapping flight,” Science, vol. 324, pp. 252– Mellon University, Pittsburgh, PA, in 2008 and 2011,
255, 2009. respectively.
[13] F. V. Breugel, W. Regan, and H. Lipson, “From insects to machines: He is currently a Research Fellow at Children’s
Demonstration of a passively stable, untethered flapping-hovering micro- Hospital Boston, Boston, MA, where he is involved
air vehicle,” IEEE Robot. Autom. Mag., vol. 15, no. 4, pp. 68–74, Dec. in the field of surgical robotics. His research interests
2008. include flapping wing flight, small-scale design and
[14] R. J. Wood, “The first takeoff of a biologically inspired at-scale robotic manufacturing techniques, and bioinspired robotics.
insect,” IEEE Trans. Robot., vol. 24, no. 2, pp. 341–347, Apr. 2008. Dr. Arabagi placed first in the 2011 ASME Grad-
[15] E. Steltz and R. Fearing, “Dynamometer power output measurements uate Robot Design Competition for his work on the design of a flapping wing
of miniature piezoelectric actuators,” IEEE/ASME Trans. Mechatronics, aerial robot.
vol. 14, no. 1, pp. 1–10, Feb. 2009.
[16] S. N. Mahmoodi, N. Jalili, and M. Daqaq, “Modeling, nonlinear dynamics,
and identification of a piezoelectrically actuated microcantilever sensor,”
IEEE/ASME Trans. Mechatronics, vol. 13, no. 1, pp. 58–65, Feb. 2008. Lindsey Hines received the B.Sc. degree in mechan-
[17] R. J. Wood, E. Steltz, and R. S. Fearing, “Optimal energy density piezo- ical engineering and the B.A. degree in mathematics
electric bending actuators,” Sens. Actuators A, vol. 119, pp. 476–488, from the University of St. Thomas, Saint Paul, MN, in
2005. 2008, and the M.Sc in robotics from Carnegie Mellon
[18] M. Sitti, D. Campolo, J. Yan, R. S. Fearing, T. Su, D. Taylor, and T. Sands, University, Pittsburgh, PA, in 2011, where she is cur-
“Development of PZT/PZN-PT unimorph actuators for micromechani- rently working toward the Ph.D. degree in robotics.
cal flapping mechanisms,” in Proc. IEEE Robot. Autom. Conf., 2001, Her interests include flapping flight and robust
pp. 3839–3846. control.
[19] M. Sitti, “Piezoelectrically actuated four-bar mechanism with two flexi- Ms. Hines has been awarded both the National Sci-
ble links for micromechanical flying insect thorax,” IEEE/ASME Trans. ence Foundation and National Defense Science and
Mechatronics, vol. 8, no. 1, pp. 26–36, Mar. 2003. Engineering Graduate Fellowship, and placed first in
[20] R. J. Wood, “Liftoff of a 60mg flapping-wing mav,” in Proc. IEEE/RSJ the 2011 Graduate Robot Design Competition.
Int. Conf. Intell. Robots Syst., 2007, pp. 1889–1894.
[21] M. Karpelson, G. Wei, and R. J. Wood, “Milligram-scale high-voltage
power electronics for piezoelectric microrobots,” in Proc. Int. Conf. Robot.
Autom., 2009, pp. 2217–2224. Metin Sitti (S’94–A’99–M’99–SM’08) received the
[22] R. J. Wood, E. Steltz, and R. Fearing, “Nonlinear performance limits for B.Sc. and M.Sc. degrees in electrical and electron-
high energy density piezoelectric bending actuators,” in Proc. IEEE/RSJ ics engineering from Bogazici University, Istanbul,
Int. Conf. Intell. Robots Syst., 2005, pp. 3633–3640. Turkey, in 1992 and 1994, respectively, and the Ph.D.
[23] S. Avadhanula, “The design and fabrication of the MFI thorax based on op- degree in electrical engineering from The University
timal dynamics,” Ph.D. dissertation, Dept. Mech. Eng., Univ. California, of Tokyo, Tokyo, Japan, in 1999.
Berkeley, CA, 2005. He was a Research Scientist at the University of
[24] L. I. Sedov, Two-Dimensional Problems in Hydrodynamics and Aerody- California, Berkeley, from 1999 to 2002. He is cur-
namics. New York: Interscience, 1965. rently a Professor in the Department of Mechanical
[25] M. H. Dickinson, F. Lehman, and S. P. Sane, “Wing rotation and the aero- Engineering and Robotics Institute, Carnegie Mellon
dynamic basis of insect flight,” Science, vol. 284, pp. 1954–1600, 1999. University, Pittsburgh, PA, where he is also the Di-
[26] A. Andersen, U. Pesavento, and Z. J. Wang, “Unsteady aerodynamic of rector of the NanoRobotics Laboratory and the Center for Bio-Robotics. He
fluttering and tumbling plates,” J. Fluid Mech., vol. 541, pp. 65–90, 2005. was the Adamson Career Faculty Fellow during 2007–2010. His research in-
[27] J. P. Whitney and R. J. Wood, “Aeromechanics of passive rotation in terests include micro/nanorobotics, bioinspired miniature mobile robots, and
flapping flight,” J. Fluid Mech., vol. 660, pp. 197–220, 2010. micro/nanomanipulation.
[28] W. B. Dickson, A. D. Straw, C. Poelma, and M. H. Dickinson, “An Dr. Sitti received the Society of Optics and Photonics Nanoengineering Pi-
integrative model of insect flight control,” in Proc. AIAA Aerosp. Sci. oneer Award in 2011. He received a National Science Foundation CAREER
Meet. Exhib., 2006, pp. 1–19. Award in 2005. From 2008 to 2010, he was the Vice President of Technical Ac-
[29] R. J. Wood, “Design, fabrication, and analysis of a 3 DOF, 3 cm flapping- tivities of the IEEE Nanotechnology Council. He was elected as a Distinguished
wing mav,” in Proc. 2007 IEEE/RSJ Int. Conf. Intell. Robots Syst., Lecturer of the IEEE Robotics and Automation Society from 2006 to 2008. He
pp. 1576–1581. also received the Best Paper Award at the IEEE/RSJ International Conference on
[30] R. Dudley, The Biomechanics of Insect Flight: Form, Function and Evo- Intelligent Robots and Systems in 2009 and 1998, the Best Biomimetics Paper
lution. Princeton, NJ: Princeton Univ. Press, 1999. Award at the IEEE Robotics and Biomimetics Conference in 2004, and the Best
[31] S. Montambault and C. M. Gosselin, “Analysis of underactuated mechan- Video Award at the IEEE Robotics and Automation Conference in 2002. He is
ical grippers,” J. Mech. Design, vol. 123, pp. 367–375, 2001. Co-Editor-in-Chief of the Journal of Micro/Nano-Mechatronics.
[32] R. Olfati-Saber, “Nonlinear control of underactuated mechanical systems
with application to robotics and aerospace vehicles,” Ph.D. dissertation,
Dept. Electr. Eng. Comput. Sci., Massachussets Inst. Technol., Cambridge,
MA, 2000.