Sage Journals: Discover world-class research

Abstract

Hovering flapping wing flight is intrinsically unstable in most cases and requires active flight stabilization mechanisms. This paper explores the passive stability enhancement with the addition of top and bottom sails, and the capability to predict the stability from a very simple model decoupling the roll and pitch axes. The various parameters involved in the dynamical model are evaluated from experiments. One of the findings is that the damping coefficient of a bottom sail (located in the flow induced by the flapping wings) is significantly larger than that of a top sail. Flight experiments have been conducted on a flapping wing robot of the size of a hummingbird with sails of various sizes and the observations regarding the flight stability correlate quite well with the predictions of the dynamical model. Twelve out of 13 flight experiments are in agreement with stability predictions.

Keywords

Micro air vehicles passive stabilization aero-dampers stability derivatives

Introduction

The complex unsteady aerodynamic mechanisms generated by insects and humming birds in hovering flight have been gradually understood over the past decades.^1–4 More recently, the extreme miniaturization of avionics stimulated the engineering community to consider building robots mimicking the behavior of insects and birds, leading to impressive projects such as Delfly,⁵ Harvard’s Robobee,⁶ Festo’s robotic Seagull,⁷ AeroVironment’s Nano Hummingbird,⁸ University of Texas A&M,⁹ or Konkuk University in Korea,¹⁰ to quote only a few. Beyond the mere curiosity of mimicking nature, it is believed that the ornithopters will one day outperform in agility the best quadcopters.

Our own project (Figure 1), named COLIBRI^11,12 flew for the first time in June 2016.¹³ In this particular design, the wings have only a single degree of freedom (flapping) and the wing shape (camber and angle of attack) is obtained passively as a result of the aerodynamic forces exerted on the wing during flapping. The wing consists of a stiffened membrane attached to two bars, the leading edge bar used for flapping and the root-edge bar which controls the aerodynamic profile of the wing during flapping. The attitude control moments are obtained by moving the root-edge bars in such a way to create a dissymmetry in the lift force distribution produced by the wings and moving the center of pressure along the wing span; a dissymmetry between the left and right wing will produce a roll moment, and a dissymmetry between the front and back half strokes will produce a pitch moment; this mechanism is known as wing twist modulation.⁸

Figure 1.

General view of the COLIBRI robot.

Most of the study reported in this paper was done before the first actively stabilized flight of the COLIBRI robot, at a time when the wing design did not generate enough lift to include all the hardware necessary for active control. The purpose of the study is to improve the understanding of the vehicle dynamics and stability in flight, which is helpful for designing the controller and achieving the actively stable flight. It is focused on the stability enhancement by means of sails, the capability to predict the stability with simple, uncoupled equations considering the robot as a rigid body, and the experimental determination of the system parameters.

Rigid body dynamics near hovering

The dynamics of a flapping wing robot near hovering may be described approximately as a rigid body; besides, the longitudinal and lateral dynamics are only weakly coupled, so that it may be assumed that they are uncoupled; they can be described by linearized Newton–Euler equations; a similar approach has been followed by e.g. Van Breugel et al.,¹⁴ Ristroph et al.,¹⁵ and Teoh et al.¹⁶ For the sake of simplicity, we will focus on the longitudinal dynamics; the lateral dynamics is similar with appropriate changes in the numerical values.

Consider the force diagram of Figure 2. At hovering equilibrium, the lift balances the weight, L = mg, and the robot is upright (pitch angle θ = 0). If a disturbance induces some θ, the thrust vector L rotates, generating a horizontal component $L \sin θ ≃ L θ$ which induces some horizontal motion (velocity u), in turn generating some opposing drag force. The horizontal velocity u modifies the wing velocity distribution w, increasing it to w + u during the upstroke and decreasing it to w − u during the downstroke (according to the coordinate system of Figure 2). If the drag force varies quadratically with the wing tip absolute velocity, the total drag force reads

f_{d} = - β {(w + u)}^{2} + β {(w - u)}^{2} = - 4 β wu = - Ku

(1)

where β is an aerodynamic constant depending on the wing shape; the damping force is linear in u. According to equation (1), the damping constant K is a linear function of the wing velocity w, that is of the flapping frequency. K may be determined experimentally with a pendulum experiment conducted with and without flapping the wing; a sketch of the experimental setup is shown in Figure 3(a) (the same set-up can be used to determine the damping constant in the lateral direction by rotating the robot by

90 °

). Figure 3(b) shows the value of K in the longitudinal and in the lateral directions for various flapping frequencies; one sees that K varies nearly linearly with the flapping frequency, as suggested from equation (1). Figure 3(c) shows typical time histories of the free response of the pendulum, with the exponential fit in dashed lines. The data presented in Figure 3 confirms the linear dependency of K on the flapping frequency; K also depends on the wing shape as illustrated in Table 1, which gives the numerical values of the longitudinal and lateral damping coefficient K of two among the many wings used in this project (the reference number in the first column is internal to the project).

Figure 2.

Coordinate system and force diagram of forces for the longitudinal (pitch) equilibrium.

Figure 3.

(a) Pendulum experiment for the determination of the damping constant $K = - X_{u}$ . (b) Damping coefficient K in the longitudinal and lateral directions for various flapping frequencies (wing MLF72E-5); the dotted lines shows the linear fit passing through the origin. (c) Typical time histories of the free response of the pendulum, with the exponential fit.

Table 1.

Robot damping coefficient K (mN.s/m) for two wings used in this project.

	Longitudinal	Lateral
MLF72-2 at 25 Hz	13.8	16.1
MLF72E-5 at 19.5 Hz	22.3	26.8

The damping forces associated with the wing motion constitute the dominant damping mechanism in the system and we will assume that all the damping forces can be reduced to a point force acting at the center of drag D and proportional to the linear velocity of the center of drag

f_{d} = - K (u + q z_{d})

(2)

where u is the velocity of the center of mass C,

q = \dot{θ}

is the pitch rate, and z_d is the distance between the center of mass and the center of drag (

z_{d} > 0

if D is above C). From Figure 2, one sees that if the center of drag is above the center of mass, the drag force generates a pitch moment which tends to reduce the pitch angle θ.

The longitudinal (pitch) and lateral (roll) dynamics may be modeled in the same way; in the following, we limit the presentation to the pitch dynamics, using classical notations of aircraft dynamics.¹⁷

Near hovering, the longitudinal dynamics is governed by Newton’s equation

m \dot{u} = X_{u} u + X_{q} q + mg θ

(3)

where the three terms in the right-hand side are respectively the drag force due to the axial velocity of the center of mass u, the drag force due to the pitch rate

q = \dot{θ}

, and the horizontal component of the wing thrust vector (assuming that θ is small and that the vertical component of the thrust equilibrates the weight mg); from equation (2),

X_{u} = - K

and

X_{q} = - K z_{d}

. The rotational equilibrium (Euler equation) reads

I \dot{q} = M_{u} u + M_{q} q + τ

(4)

where I is the moment of inertia about the center of mass C and the three terms on the right-hand side are respectively the moment with respect to C of the drag forces associated with the translation velocity u, the moment of the drag forces due to the rotation velocity q, and the external control torque, τ. According to our assumption that the drag forces can be reduced to a point force acting at the center of drag,

M_{u} = - K z_{d} = X_{q}

. However, for the rotational damping, we include an additional term K_r representing the rotational damping about the center of drag, so that

M_{q} = - K_{r} - K z_{d}^{2}

. M_q can be estimated with a pendulum experiment similar to that used to evaluate the translational damping, with the bird attached to the center of mass (Figure 4(a)). Figure 4(b) shows the magnitude of M_q measured for different flapping frequencies; M_q varies nearly linearly with the flapping frequency. Figure 4(c) shows typical time histories of the free response of the pendulum, with the exponential fit in dashed lines.

Figure 4.

(a) Pendulum experiment for the determination of the rotational damping constant M_q. (b) Magnitude of the damping coefficient $M_{q} (M_{q} < 0)$ appearing in the longitudinal equilibrium equation for various flapping frequencies (wing MLF72E-5); the dotted lines shows the linear fit passing through the origin. (c) Typical time histories of the free response of the pendulum, with the exponential fit.

The equations may be written in state space, using the state vector ${(u, q, θ)}^{T}$

{\begin{matrix} \dot{u} \\ \dot{q} \\ \dot{θ} \end{matrix}} = [\begin{matrix} {\hat{X}}_{u} & {\hat{X}}_{q} & g \\ {\hat{M}}_{u} & {\hat{M}}_{q} & 0 \\ 0 & 1 & 0 \end{matrix}] {\begin{matrix} u \\ q \\ θ \end{matrix}} + {\begin{matrix} 0 \\ 1 \\ 0 \end{matrix}} \frac{τ}{I}

(5)

where

{\hat{X}}_{u} = X_{u} / m

and

{\hat{M}}_{q} = M_{q} / I

are always negative and

{\hat{X}}_{q} = - K z_{d} / m

and

{\hat{M}}_{u} = - K z_{d} / I

are negative if

z_{d} > 0

, that is if the center of drag is above the center of mass, and positive if

z_{d} < 0

The characteristic equation reads

λ^{3} - ({\hat{X}}_{u} + {\hat{M}}_{q}) λ^{2} + ({\hat{M}}_{q} {\hat{X}}_{u} - {\hat{M}}_{u} {\hat{X}}_{q}) λ - {\hat{M}}_{u} g = 0

(6)

The pole pattern obtained with this fairly simple model is consistent with more elaborate models available in the literature.^18–22

Passive stability with sails

Consider the system of Figure 5 where the flapping wing robot has been supplemented with two sails, a top sail of area S₁ at z₁ above the center of mass and a bottom sail of area S₂ at z₂ below the center of mass. Again, we assume that the drag forces acting on the sails can be reduced to point forces acting at the geometrical center of the sail, and proportional to the absolute linear velocity of the geometrical center: $f_{i} = - k_{i} S_{i} v_{i}$ ; this linear viscous damping assumption is confirmed by the experiments reported in the next section. The whole system is once again considered as a rigid body of mass m and center of mass C. With these assumptions, the various terms involved in the longitudinal dynamics are as follows:

Figure 5.

Robot configuration with top and bottom sails (z₁ and z₂ are both positive; z_d is positive if the center of drag D is above the center of mass C).

Drag force due to the axial velocity u

X_{u} u = - Ku - k_{1} S_{1} u - k_{2} S_{2} u = (- K - k_{1} S_{1} - k_{2} S_{2}) u

(7)

Pitch moment due to u

M_{u} u = (- K z_{d} - k_{1} S_{1} z_{1} + k_{2} S_{2} z_{2}) u

(8)

Drag force due to the rotational velocity q

X_{q} q = (- K z_{d} - k_{1} S_{1} z_{1} + k_{2} S_{2} z_{2}) q

(9)

Pitch moment due to the rotational velocity q

M_{q} q = (- K z_{d}^{2} - k_{1} S_{1} z_{1}^{2} - k_{2} S_{2} z_{2}^{2}) q

(10)

Notice that, once again, M_u = X_q. We will address shortly how the drag coefficients of the sails may be determined experimentally. Before this, let us consider the inertia properties of the system. The vertical equilibrium of the system implies that $L = mg \cos θ ≃ mg$ while the horizontal component of the wing thrust force is $mg \sin θ ≃ mg θ$ .

Let m₁ be the added mass of air associated with the top sail; the mass of air which can be regarded as the mass attached to the sail (e.g. see White²³), and m₂ the added mass of air of the bottom sail, that we assume lumped at the center of the sails, respectively at z₁ above C and z₂ below C. Newton’s equation describing the longitudinal dynamic equilibrium becomes

(m + m_{1} + m_{2}) \dot{u} = X_{u} u + X_{q} q + mg θ

(11)

and Euler’s equation which describes the pitch equilibrium about θ becomes

(I + m_{1} z_{1}^{2} + m_{2} z_{2}^{2}) \dot{q} = M_{u} u + M_{q} q

(12)

where I is the total moment of inertia about C and the various terms involved in the right-hand side are defined by equations (7)–(10). Notice that the added masses appear only in the inertia forces and not in the horizontal component of the thrust force. The foregoing equations may be casted in a state-space form similar to equation (5)

{\begin{matrix} \dot{u} \\ \dot{q} \\ \dot{θ} \end{matrix}} = [\begin{matrix} {\hat{X}}_{u} & {\hat{X}}_{q} & g^{*} \\ {\hat{M}}_{u} & {\hat{M}}_{q} & 0 \\ 0 & 1 & 0 \end{matrix}] {\begin{matrix} u \\ q \\ θ \end{matrix}}

(13)

where

{\hat{X}}_{u} = X_{u} / (m + m_{1} + m_{2})

and

{\hat{M}}_{q} = M_{q} / (I + m_{1} z_{1}^{2} + m_{2} z_{2}^{2})

are always negative,

{\hat{X}}_{q} = X_{q} / (m + m_{1} + m_{2})

and

{\hat{M}}_{u} = M_{u} / (I + m_{1} z_{1}^{2} + m_{2} z_{2}^{2})

, and

g^{*} = mg / (m + m_{1} + m_{2})

. Since the sails have the same shape in pitch and roll, the equations are very similar for the lateral (roll) axis, except for a different value of the damping constant K (Table 1).

Added mass and damping of the sails

The sails appear in the foregoing equations through the viscous damping and the added mass (the mass of air which is moving with the sail), respectively $k_{1} S_{1}$ and m₁ for the top sail and $k_{2} S_{2}$ and m₂ for the bottom sail. Surprisingly, it has been observed that the behavior of the bottom sail is significantly different from that of the top sail, because of the downflow induced by the flapping wings.

To evaluate the added mass and the damping coefficient, a vibration experiment was conducted in which the sail is attached to a cantilever beam (Figure 6); the beam is excited by a voice coil and the beam vibration is monitored with a laser vibrometer; the frequency response functions (FRFs) are recorded, first without sails and then with sails of various sizes. Additionally, a flapping wing mechanism is used to simulate the downward flow when studying the bottom sail. A finite element model of the system (cantilever beam + point mass and damper at the geometrical center of the sail) has been developed and the damping coefficient and the added mass of every sail are calculated by curve fitting on the FRFs (Figure 7). The good quality of the fit confirms the assumptions made.

Figure 6.

Experimental setup for the determination of the added mass and damping of the sails. The sail is attached to the tip of a cantilever beam. A flapping wing mechanism (not used for the top sail) is mounted on a separate support to simulate the air flow acting on the bottom sail. d is the distance between the geometrical center of the sail and the wing root.

Figure 7.

FRF of the system of Figure 6 for a top sail of various sizes (full line: experiment; dotted line: finite element model including an added point mass and damper at the geometric center of the sail).

The damping coefficient $k_{1} S_{1}$ and the added mass of air m₁ of the top sail are reported in Table 2. For the bottom sail, the flapping wing mechanism is used at the normal flapping frequency (21 Hz) at a distance d between the geometrical center of the sail and the wing root (corresponding to the configurations used in the flight experiments). The experiment led to the surprising results of Figure 8 that the added mass is not significantly affected by the air flow while the damping coefficient is. Due to the down flow induced by the flapping wings, the damping coefficient measured for the bottom sail is one order of magnitude larger than that of the top sail and depends critically on the distance d between the geometrical center of the sail and the wing root. The damping coefficient $k_{2} S_{2}$ of a bottom sail of 50 cm² is reported in Table 3 for two values of the distance d. The above data allow to compute all the parameters of the linearized model of the foregoing section.

Table 2.

Top sail damping coefficients and added mass of air.

Sail surface (cm²)	$k_{1} S_{1}$ ( $mN . s / m$ )	m₁ (g)
50	1	0.3
150	5	1.1
220	5	2.8
600	38	10.0

Figure 8.

FRF of the system of Figure 6 for a bottom sail of 50 cm² at various distances from the wing root ( $d =$ 10 cm and 23 cm); wing MLF72-2 flapping at 21 Hz (full line: experiment; dotted line: finite element model).

Table 3.

Bottom sail of 50 cm²: damping coefficient and added mass of air as a function of the distance d between the geometrical center of the sail and the wing root, for wing MLF72-2 flapping at 21 Hz.

d (cm)	$k_{2} S_{2}$ ( $mN . s / m$ )	m₂ (g)
10	15	0.29
23	18	0.32

Flight experiments

Thirteen flight experiments have been conducted with the flapping wing robot equipped with top and bottom sails of various sizes (Figure 9); for every configuration, the position of the center of mass C was determined experimentally; the position of the center of drag D with respect to the center of mass is calculated by the formula $z_{d}^{*} = M_{u} / X_{u}$ .

Figure 9.

Flapping wing robot equipped with top and bottom sail stabilizers used in the flight experiments. The reflective markers used for motion tracking in the video room are also indicated.

Figure 10 and Table 4 describe the various configurations: mass, size of the sails, $z_{d}^{*}$ . Table 5 gives the numerical values of the stability derivatives in pitch [components of the system matrix, equation (13)] and in roll. Table 6 gives the eigenvalues in pitch and roll and the predicted behavior: I-D means “Instable-Divergent” (one positive real eigenvalue); I-O means “Instable-Oscillatory” (a pair of complex eigenvalues with positive real part); S means “Stable.” The last column of the table gives the behavior observed during the flight; videos showing the various behaviors observed are shown in Video of flight with sails.²⁴

Figure 10.

Configurations used in the flight experiments. The weight of the robot equipped with sails in each configuration is shown below it.

Table 4.

Characteristics of the various flight configurations: mass, top and bottom sail area and position, added mass of air, position of the global center of drag $z_{d}^{*}$ .

Flight N ^o	M (g)	S₁ (cm²)	m₁ (g)	z₁ (cm)	S₂ (cm²)	m₂ (g)	z₂ (cm)	$z_{d}^{*}$ (cm)
1	16.3	600	10	18.3	–	–	–	13
2	17.0	600	10	15.2	–	–	–	10
3	17.0	600	10	19.8	–	–	–	14
4	17.1	600	10	19.0	50	0.1	12.8	8
5	17.5	600	10	19.8	50	0.3	25.0	4
6	18.2	600	10	18.3	50	0.3	26.5	3
7	17.0	600	10	11.2	–	–	–	6
8	16.7	600	10	7.50	–	–	–	4
9	16.8	150	1.1	20.2	50	0.1	14.7	−6
10	16.8	150	1.1	23.9	50	0.1	11.0	−2
11	17.2	150	1.1	24.8	50	0.2	22.6	−8
12	16.9	150	1.1	20	50	0.2	17.7	−4
13	17.0	200	1.8	15.5	50	0.1	10.0	0

Table 5.

Stability derivatives of longitudinal (column 2–5) and lateral (column 6–9) dynamics for 13 flight tests with the configurations shown in Figure 10.

Flight N ^o	${\hat{X}}_{u}$	${\hat{X}}_{q}$	${\hat{M}}_{u}$	${\hat{M}}_{q}$	${\hat{Y}}_{v}$	${\hat{Y}}_{p}$	${\hat{L}}_{v}$	${\hat{L}}_{p}$	$g *$
1	−2.05	−0.26	−15.56	−3.09	−2.13	−0.26	−15.37	−3.09	6.08
2	−1.99	−0.19	−13.32	−2.54	−2.08	−0.19	−12.93	−2.56	6.18
3	−1.99	−0.28	−15.09	−3.11	−2.08	−0.28	−15.00	−3.11	6.18
4	−2.54	−0.20	−11.07	−3.56	−2.62	−0.19	−10.94	−3.56	6.19
5	−2.61	−0.11	−5.71	−4.90	−2.69	−0.11	−5.63	−4.90	6.24
6	−2.55	−0.07	−3.77	−4.77	−2.63	−0.07	−3.62	−4.78	6.33
7	−1.99	−0.11	−7.02	−1.53	−2.08	−0.10	−6.46	−1.59	6.18
8	−2.02	−0.09	−10.21	−1.15	−2.10	−0.08	−9.71	−1.17	6.14
9	−1.89	0.11	8.61	−2.54	−2.02	0.11	9.15	−2.57	9.21
10	−1.89	0.04	4.38	−3.11	−2.02	0.04	4.61	−3.12	9.21
11	−2.01	0.16	15.79	−6.66	−2.14	0.16	15.86	−6.66	9.22
12	−2.05	0.09	10.36	−5.14	−2.17	0.08	9.72	−5.17	9.21
13	−1.85	−0.01	−1.29	−3.45	−1.97	−0.01	−1.17	−3.45	8.87

Table 6.

Eigenvalues of longitudinal (column 2 and 3) and lateral (column 4 and 5) dynamics versus the observations.

Flight N ^o	λ ₁	$λ_{2} and λ_{3}$	λ ₄	$λ_{5} and λ_{6}$	Predicted stability	Observation
1	−6.83	$0.85 \pm 3.62 j$	−6.84	$0.81 \pm 3.61 j$	I–O	I–O
2	−6.24	$0.85 \pm 3.53 j$	−6.24	$0.80 \pm 3.49 j$	I–O	I–O
3	−6.82	$0.86 \pm 3.60 j$	−6.84	$0.82 \pm 3.59 j$	I–O	I–O
4	−6.63	$0.26 \pm 3.20 j$	−6.63	$0.22 \pm 3.19 j$	I–O	S
5	−6.48	$- 0.51 \pm 2.29 j$	−6.49	$- 0.55 \pm 2.26 j$	S	S
6	−6.00	$- 0.66 \pm 1.88 j$	−5.99	$- 0.71 \pm 1.82 j$	S	S
7	−4.87	$0.68 \pm 2.90 j$	−4.83	$0.58 \pm 2.82 j$	I–O	I–O
8	−5.21	$1.02 \pm 3.32 j$	−5.17	$0.95 \pm 3.26 j$	I–O	I–O
9	3.02	$- 3.72 \pm 3.53 j$	3.06	$- 3.82 \pm 3.60 j$	I–D	I–D
10	2.02	$- 3.51 \pm 2.75 j$	2.04	$- 3.60 \pm 2.81 j$	I–D	I–D
11	3.09	$- 5.88 \pm 3.55 j$	3.05	$- 5.93 \pm 3.57 j$	I–D	I–D
12	2.66	$- 4.93 \pm 3.40 j$	2.53	$- 4.93 \pm 3.32 j$	I–D	I–D
13	−4.45	$- 0.43 \pm 1.55 j$	−4.42	$- 0.50 \pm 1.45 j$	S	S

Figure 11 shows examples of pitch and roll signals recorded during flight tests. In order to capture the attitude, a room equipped with 8 OptiTrack Flex motion tracking cameras has been used. These cameras offer the resolution of 1280 px × 1024 px and the frequency of 120 frames per second (FPS). Five retroreflective markers were attached to the robot with a dissymmetric arrangement to avoid the loss of orientation and attitude from the tracking system. In analyzing Table 6, it is interesting to note that all the observed behaviors during the flights are in agreement with the predictions except for flight $N^{o} 4$ . Figure 12 shows the real part of the eigenvalues of the longitudinal dynamics predicted by our model, as a function of the distance between the center of mass and the center of drag, $z_{d}^{*}$ ( $z_{d}^{*}$ is positive when the center of drag is above the center of mass). We note that:

Figure 11.

Examples of pitch (left) and roll (right) signals recorded during flight tests in the video tracking room. From top to bottom: Stable flight ( $N^{o} 13$ ), oscillatory unstable ( $N^{o} 3$ ), and divergent ( $N^{o} 11$ ).

Figure 12.

Numerical predication of the real part $(σ)$ of the eigenvalues of the longitudinal dynamics as a function of $z_{d}^{*}$ . All flight tests confirm the numerical predictions except for flight $N^{o} 4$ which was predicated unstable and observed stable.

None of the flights with top sail only was stable.

All the twin sails flights with negative $z_{d}^{*}$ were instable.

All stable cases had a small positive value of $z_{d}^{*}$ .

All predicted stable flights were observed stable with the exception of flight $N^{o} 4$ , which was very close to the stability limit.

The same conclusions apply to the lateral dynamics.

Conclusion

This paper has analyzed the dynamic stability of a flapping twin-wing robot near hovering; a very simple model (similar to those used in aircraft dynamics) decoupling pitch and roll has been used to show that the system is intrinsically unstable. The model has been used to study the passive stability enhancement with the addition of top and bottom sails. Experiments have been conducted to evaluate the parameters involved in the dynamical model; the experiments revealed that the damping coefficient of the bottom sail (located in the flow induced by the flapping wings) is significantly larger than that of the top sail and depends critically on the flapping frequency and the distance between the geometrical center of the sail and the wing root. Thirteen flight experiments have been conducted with sails of various sizes and the behavior of the robot was observed; 12 out of 13 flight experiments are in agreement with stability predictions of our simplified model. This led to trust the model and use it later for designing the controller, and achieving the actively stable flight. The study indicates that $z_{d}^{*}$ plays an important role on stability; none of the flights with negative values of $z_{d}^{*}$ , nor with large positive values (typical of single sail configuration) were stable. In spite of the variety of sail sizes, the model was able to successfully predict the stability.

Footnotes

Acknowledgements

The authors wish to thank the contribution of Ignacio Senet Capote, Mylène Dumon, Carlos Santos, and Georges Tod.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

References

Weis-Fogh

Unusual mechanisms for the generation of lift in flying animals.

Scient Am 1975; 233: 80–87.

Ellington

CP.

The novel aerodynamics of insect flight: application to micro-air vehicles. J Exp Biol 1999; 202: 3439–3448.

Dickinson

Lehman

F-O

Sane

SP.

Wing rotation and the aerodynamic basis of insect flight.

Science 1999; 284: 1954–1960.

Sane

SP.

The aerodynamic of insect flight.

J Exp Biol 2003; 206: 4191–4208.

Lentink

Jongerius

Bradshaw

NL.

The scalable design of flapping micro-air vehicles inspired by insect flight, Ch.14. In: Floreano

et al . (eds) Flying insects and robots. Berlin Heidelberg: Springer-Verlag, 2009.

Wood

RJ.

The first takeoff of a biologically inspired at-scale robotic insect.

IEEE Trans Robot 2008; 24: 341–347.

Festo Corporation, www.festo.com/group/en/cms/10238.htm

Keennon

Klingebiel

Won

et al . Development of the Nano hummingbird: a tailless flapping wing micro air vehicle. AIAA paper 2012-0588, 2012, pp. 1–24.

Coleman

Benedict

Hrishikeshavan

et al . Design, development and flight-testing of a robotic hummingbird, AHS 71st Annual Forum, Virginia Beach, 5–7 May 2015.

10.

Phan HV, Kang T and Park HC. Design and stable flight of a 21 g insect-like tailless flapping wing micro air vehicle with angular rates feedback control. Bioinspiration & biomimetics 2017; 12(3): 036006.

11.

Karásek

Robotic hummingbird: design of a control mechanism for a hovering flapping wing micro air vehicle. PhD Thesis, Université Libre de Bruxelles, Active Structures Laboratory, Belgium, 2014.

12.

Roshanbin

Altartouri

Karásek

et al . COLIBRI: a hovering flapping twin-wing robot. Int J Micro Air Vehicles 2017; 9: 270–282.

13.

Video of actively stabilized flight (Click here to visualize: https://www.youtube.com/watch?v=aWeUPiz2pt4).

14.

Van Breugel

Regan

Lipson

From insects to machines.

IEEE Robot Automat Mag 2008; 15: 68–74.

15.

Ristroph

Morozova

et al . Active and passive stabilization of body pitch in insect flight. JR Soc Interface 2013; 10: 20130237.

16.

Teoh

Fuller

Chirarattananon

et al . A hovering flapping-wing microrobot with altitude control and passive upright stability. In: IEEE/RSJ international conference on intelligent robots and systems, Vilamoura, 2012, pp. 3209–3216.

17.

Franklin

JA.

Dynamics, control and flying qualities of V/STOL aircraft. Reston, VA: AIAA Education Series, 2002.

18.

Sun

Xiong

Dynamic flight stability of a hovering bumblebee.

J Exp Biol 2005; 208: 447–459.

19.

Orlowski

Girard

AR.

Dynamics, stability, and control analyses of flapping wing micro-air vehicles.

Prog Aerospace Sci 2012; 51: 18–30.

20.

Zhang

Sun

Dynamic flight stability of a hovering model insect: lateral motion.

Acta Mech Sin 2010; 26: 175–190.

21.

Karásek M and Preumont A. Flapping flight stability in hover: A comparison of various aerodynamic models. International Journal of Micro Air Vehicles 2012; 4(3): 203–226.

22.

Taha

Hajj

Nayfeh

AH.

Longitudinal flight dynamics of hovering MAVs/insects.

J Guid Contr Dyn 2014; 37: 970–978.

23.

White

FM.

Fluid mechanics. New York, NY: McGraw-Hill, 7th edition, 2011.

24.

Video of flight with sails (Click here to visualize: https://www.youtube.com/watch?v=w6tFYjUCSIM).