Effect of wing flexibility on flight dynamics stability of flapping wing MAVs in forward flight

Abstract

The paper presents a qualitative study of the effect of the wing flexibility on the flight stability of flapping wing micro air vehicles in forward flight. The longitudinal dynamic flight stability is compared between rigid wing flapping wing micro air vehicles and flexible-wing flapping wing micro air vehicles. The aerodynamics derivatives are computed respectively using the method of computational fluid dynamics and fluid–structure interaction method, and the techniques of eigenvalue and eigenvector analysis are applied to solve the equations of motion. It is shown that the flexibility can change the stability properties of the flapping micro air vehicles from unstable to stable. Besides, the position of mass center of flapping micro air vehicles is also a key factor to the flight stability, which provides a new insight to design a stable flapping wing micro air vehicle.

Introduction

In the last 20 years, flapping wing micro air vehicle (FMAV) has aroused intense interests from researchers because of their potential to dramatically strengthen our sensing and information gathering capabilities. ¹ In order to design a series of FMAVs with high aerodynamic efficiency and maneuverability, researchers have devoted much effort to handle the challenging problems, including unsteady aerodynamics, aero-elasticity, nonlinear flight dynamics, etc. Several FMAVs have already been designed and fabricated successfully to accomplish the military and civilian mission, such as the American Hummingbird and Chinese Dove. ²

Insects and birds in nature with outstanding flying skills in a broad range of low Reynolds number provide best inspiration for researchers to develop the FMAV. Taylor and Thomas ³ made the first formal quantitative research of dynamic stability of the desert locust using the linear six-DOF dynamic model, whose aerodynamic derivatives obtained experimentally. However, this experimental approach had to depend on some control responses by insects. In order to avoid such problem, Sun and Xiong ⁴ used the computational fluid dynamics (CFD) method to compute the aerodynamics derivatives, which were applied for eigenvalue and eigenvector analysis to study the longitudinal and lateral stability in several hovering insects.^4,5

In the above analysis, it is notable that both studies are based on two assumptions. One is that the flapping wing system was assumed to be a time-invariant model through neglecting the wing mass and averaging the aerodynamic forces over a wing flapping cycle (averaged model). In order to consider the effect of the periodic effect and nonlinearity of the aerodynamic forces, a nonlinear rigid body model of an insect was developed by Bierling and Patil. ⁶ The stability of the linearized, time-averaged system was first calculated and then the stability of the time-variant system was studied using Floquet theory. When the wing mass was less than 1% of the total mass, both approaches obtained similar results. Sun et al. ⁵ obtained same conclusion when they used same method to study two hovering model of insects: a dragonfly and a hawk moth. They concluded the model with relatively high wing-beat frequency and small wing-mass to body-mass ratio is suitable for the averaged model theory. In the present paper, the research model is Chines Dove, as shown in Figure 1, whose wing-beat frequency is around 10–20 Hz and wing mass-to-body mass ratio is 2.0%. Table 1 presents the flapping frequency and forward velocity of DelFly II and Cockatoo of the same size in comparison with Chinese Dove. According to Caetano et al., ⁷ if the single wing-mass to body-mass ratio of FWMAVs is under 2.8%, application of single rigid body kinematics for the dynamic simulation can ensure that considerable information is not lost. Consequently, Chinese Dove can be approximately treated as an averaged model and linearized system.

Figure 1.

DOVE FMAV.

Table 1.

Some classical FMAV and bird species of the same size in comparison with Chinese Dove.

Parameters	DelFly II	Cockatoo	Chinese Dove
Flapping frequency	8 Hz	7.5 Hz	10–20 Hz
Averaged forward speed	7 m/s	4.4 m/s	10 m/s
Wing-to-body mass ratio	1.8%	7%	2%

Another assumption is that the wing is rigid. There is no formal quantitative study or even qualitative discussion about the effect of wing deformation on flight dynamics stability. Despite Sun ⁸ discussed the effect of wing deformation on flight dynamics for insects and concluded that relatively small deformation produces negligibly small variations in aerodynamic force during the disturbed motion and hence cannot influence the stability properties. However, by observing birds flight, it can be seen that they can actively and dramatically change the camber of wings to make the best use of aerodynamics. ⁹ At the same time, since using lightweight structures, there is also relatively large passive wing deformation for bird-like FMAV compared with the insects.^10–15 In conclusion, the influence of wing flexibility on flight stability cannot be neglected as insects.

In the present paper, the longitudinal dynamic flight stability of Chinese Dove is studied. By the method of CFD and FSI, the aerodynamic derivatives of FMAV with rigid wing and flexible wing are obtained respectively. The remainder of the paper is organized as follows: First, conditions for equilibrium of two different FMAVs in forward flight were calculated. Second, the aerodynamic derivatives of two different FMAVs at equilibrium flight were obtained. Then, employing the method of eigenvalue and eigenvector analysis, the longitudinal dynamic flight stability of two different FMAVs was analyzed. Finally, the influence of position of mass center on flight stability is also discussed and studied.

Methodology

Materials

The wing of FMAV is a quarter-elliptic platform, as shown in Figure 1. Wings are manufactured by unidirectional carbon fiber forming the wing skeleton and polyester film used as skin. The geometric and elastic parameters of each part are shown in Figure 2 and Table 2. L denotes the length of carbon rods; D denotes the diameter of rods; E is the Young’s modulus; G is the shear modulus of cross-section; ρ is the density of rods. E and G are experimental data obtained by American Instron Material Testing Equipment.

Figure 2.

Geometry of flapping wing (mm).

Table 2.

The geometric parameters and elastic properties of wing.

Part	L (mm)	D (mm)	E (GPa)	G (GPa)	ρ(g/cm3)
Fore-beam	277	2.0	120.8 ± 19.4	4.8 ± 0.1	1.55
Rear-beam	259	1.5	119.4 ± 5.1	4.4 ± 0.7	1.51
Rib	72.5–100	1.0	123.3 ± 3.4	4.0 ± 0.5	1.51
Membrane	/	/	0.15	/	1.15

Equations of motion

Similar to Sun and Xiong, ⁴ the rigid body approximation is presumed: in the present case of symmetric longitudinal motion, the FMAV is treated as a rigid body of three degrees of freedom and the cyclic forces and moments of the flapping motion are replaced by the wing-flapping-cycle average aerodynamic forces, moments and inertial forces, while assuming the gyroscopic forces of wing to be negligible.

According to Taylor and Thomas, ³ the origin o of a non-inertial body coordinate system o-x-y-z is at the center of mass of the FMAV and the pitch angle between the horizontal and the x-axis is represented by θ. The forward (u) and dorso-ventral (w) components of velocity are along x- and z- axes, respectively, and the pitching angular velocity around the center of mass is represented by q.

The inherently non-linear equations of motion could be linearized by slightly disturbing the FMAVs motion supposing that FMAV is in a steady forward flight and flight condition is symmetric about x-z plane. The linearized results ¹⁶ are:

$[\delta \stackrel{·}{u}\delta \stackrel{·}{w}\delta \stackrel{·}{q}\delta \stackrel{·}{\theta }]=[\frac{X_u}{m}\frac{X_w}{m}\frac{X_{\text{q}}}{m}-w_e-g\cos {\theta }_e\frac{Z_u}{m}\frac{Z_w}{m}\frac{Z_{\text{q}}}{m}+u_e-g\sin {\theta }_e\frac{M_u}{I_{\mathrm{yy}}}\frac{M_w}{I_{\mathrm{yy}}}\frac{M_q}{I_{\mathrm{yy}}}00010][\delta u\delta w\delta q\delta \theta ]$

(1)

where X_u, X_w, X_q, Z_u, Z_w, Z_q, M_u, M_w, and M_q are the aerodynamic derivatives about each variables. And X and Z denotes the x- and z-components of the aerodynamic force of the wing assuming the body’s aerodynamic effects negligible, respectively, and M represents the aerodynamic pitching moment; m is the mass of the FMAV; g represents the gravitational acceleration; $I_{\mathrm{yy}}$ is the pitching moment of inertia about y axis; u and w denote the x and z component of the velocity of the mass center of the FMAV, respectively; the symbol δ represents a small disturbance quantity.

In order to analyze the system conveniently, we set the pitch rate derivatives to be zero, which is discussed and proved to be feasible conditionally by Taylor and Thomas ³ Consequently, the longitudinal system matrices are only composed of the static stability derivatives (i.e.X_u, X_w, Z_u, Z_w, M_u, and M_w).

In forward flight, the equations can be further simplified because w, q, θ can be set to zero. θ is zero since the x-axis can be aligned with horizontal, and the forces and moments are in equilibrium so X = 0, Z = –mg and M = 0. With these simplifications, we may write the equation of motion for non-maneuvering flight as

$[\delta \stackrel{·}{u}\delta \stackrel{·}{w}\delta \stackrel{·}{q}\delta \stackrel{·}{\theta }]=[\frac{X_u}{m}\frac{X_w}{m}0-g\frac{Z_u}{m}\frac{Z_w}{m}{u}_{e}0\frac{M_u}{I_{\mathrm{yy}}}\frac{M_w}{I_{\mathrm{yy}}}000010][\delta u\delta w\delta q\delta \theta ]$

(2)

Flow solver

By solving the three-dimensional unsteady compressible Navier–Stokes equations, ¹⁷ the unsteady aerodynamics of the flapping wing can be derived as

$\frac{D}{\mathrm{Dt}}\int \int {\int }_{\Omega }[W]\text{d}V+\int {\int }_{\partial \Omega }\overline{\overline{H\text{'}}}·n\text{d}S=\int {\int }_{\partial \Omega }\overline{\overline{H_{\nu }}}·n\text{d}S$

(3)

where [W] is the state vector of conservative variables; Ω is the control volume; S, n are the boundary of control volume and its unit-normal outer vector, respectively; and $H\text{'}$ , H_v represent inviscid and viscous fluxes, respectively.

To develop a robust and efficient aerodynamic analysis method suitable for design of wing and complex configuration, Han et al. ¹⁸ presented and studied a combination of matrix preconditioning and multigrid method. Choi and Merkle ¹⁹ preconditioning matrix was used to precondition the time derivatives of three-dimensional Navier–Stokes equations. A cell-centered finite-volume method was used to solve the governing equations and time-stepping method utilized a lower upper symmetric-gauss-Seidel (LU-SGS) method ²⁰ with multigrid acceleration, ²¹ and turbulence model employed here is the k-ω shear stress transport (SST) model. ²²

Because there is a large scale of movement and deformation for the flexible flapping wing flight, it needs to develop an automatic mesh generation code. Based on the infinite interpolation, a moving grid methodology is applied for automatic mesh generation of the flexible flapping motion. A CO-type body-fitted grid is generated and used for the aerodynamic calculation. The outer boundary of grid is set at 20 chord lengths in the spanwise direction. Moderate and fine meshes with 0.59 $\times 106$ and 1.18 $\times 106$ grid nodes are used for grid test, respectively. To verify the reasonability of time steps, the moderate-grid simulations are performed using 72 and 144 time steps per period, respectively. Figure 3 shows lift coefficients of above aforementioned three cases when the angle of attack equals 4°, forward speed is 4 m $\text{s}-1$ , and the flapping amplitude is 70°. Through comparison, though there is little discrepancies between three cases, the moderate grid and time steps provide a sufficiently good accuracy for the simulation considering the computational cost, especially for the calculation of averaged aerodynamic forces. Therefore, the moderate spatial and temporal resolutions are used in the following studies. Figure 4 shows the grid topology of the flapping wing generated by the grid generation code.

Figure 3.

Test of grid and time steps per period.

Figure 4.

Topology of the grid distribution of the wing.

Structure solver

Based on Timoshenko beam theory, a kind of beam element is applied to simulate the wing skeleton made from carbon fiber. The global mass matrix [M], the global stiffness matrix [K], and the global force vector [F] can be obtained by assembling the element connectivity matrix, the elemental matrices, and force vectors. The equations of motion can be presented as

$[M]\{\stackrel{\mathrm{··}}{q}\}+[K]\{q\}=\{F(t)\}$

(4)

A Newmark solution is used to solve equation (4). ²³

As mentioned in Xue et al., ²⁴ there is no prestress on the membrane, which approximates that only the skeleton bears the bending loads, so it is feasible to only consider the wing skeleton during structural modeling and distribute the mass of membrane on the skeleton to consider the inertial loads of it. Table 3 shows the parameters of finite element model of wing skeleton.

Table 3.

Parameters of finite element model of wing skeleton.

FEA model	Wing skeleton
Element type	B32 (a spatial beam that uses quadratic interpolation)
Element number	130
Node member	256
Cross	Circular
Boundary condition	Fixed at two root points

If the structural damping effect is ignored, the transverse displacement y(x,t) of an elastic cantilever beam, with Young modulus E, moment of inertia I, density ρ, cross-section A, and external applied loads q(x,t), satisfies the equation of bending waves

$\frac{\partial }{\partial x2}(\mathrm{EI}\frac{\partial 2y}{\partial x2})+\rho A\frac{\partial 2y}{\partial t2}=q(x,t)$

(5)

When the beam root is prescribed as a sinusoidal angular velocity, like flapping wing, the inertial loads can be expressed as

$q(x,t)=\mathrm{kx}\cos ({\omega }_{e}t)$

(6)

Assuming transverse vibration with uniform cross-section is the linear combination of all the natural modes

$y(x,t)=\sum _{i=1}^n{Y}_i(x){\eta }_i(t)$

(7)

where ${\eta }_i(t)$ is the normal coordinate, and substituting equation (7) in equation (5) leads to

$\sum _{i=1}^n({\mathrm{EIY}}_i\text{'}\text{'})\text{'}\text{'}{\eta }_i(t)+\rho A\sum _{i=1}^n{Y}_i{\stackrel{\mathrm{··}}{\eta }}_i=q(x,t)$

(8)

By multiplying both sides of equation by Y_i, and integral x in [0,l], function can be transformed to

${\stackrel{\mathrm{··}}{\eta }}_i+{\omega }_i^2{\eta }_i=q_i(t)$

(9)

where

$q_i(t)={\int }_0^{l}q(x,t)Y_i(x)\text{d}x=h\cos ({\omega }_{\mathrm{ex}}t)$

(10)

In the case of y(x),0 = 0, $\frac{\partial y}{\partial t}(x,0)=0$ , ${\eta }_j=0$ , ${\stackrel{·}{\eta }}_j=0$ . When ${\omega }_j\ne {\omega }_e$ , the specific solution of the differential equation can be obtained

${\eta }_i=\frac{h}{{\omega }_i^2-{\omega }_e^2}(\cos ({\omega }_{e}t)\cos ({\omega }_{i}t))$

(11)

The transverse vibration can be expressed as

$y(x,t)=\sum _{i=1}^n{Y}_i(x){\eta }_i(t)$

(12)

Fluid–structure interaction

The method for fluid–structure interaction (FSI) is same as that in Yang et al. ¹⁵ Yang et al.^25–27 has focused on the flapping wing aerodynamic simulation through developing a three-dimensional unsteady Navier–Stokes solver, which is also proved to be an effective and robust tool in FMAV aerodynamic simulation, especially for the flexible wing.

The geometry is discretized in a different way for the CFD and CSD. The radial basis functions (RBFs) are applied for the information interaction between the fluid and structure solvers. The brief description of FSI process is given as follows: first the periodic aerodynamic force is obtained by fluid solver; Second, load the periodical aerodynamic force on the structure, and periodic structural deformation can be simulated by structural dynamic solver; then add the periodical deformation on the original mesh, then the periodic aerodynamic force is simulated again until both the structural deformation and the aerodynamic forces between two continuous time steps are convergent. It was tested and validated by comparing with experimental results in wind tunnel. ¹⁵ The calculated lift coefficients of flapping wing correlate well with the experimental results. The trend of computational thrust coefficients is consistent with the experimental and can be amended by taking account of the effects of wind tunnel wall interference as well as support interference on aerodynamic force. Hence, the FSI solver can be approximately considered reliable.

Force and moment equilibrium

As is known to us, the wing flexibility within a certain extent can dramatically improve the thrust force.^12,13 Thus, FMAVs with rigid wing and flexible wing correspond different cruising speed when the flapping frequency is 10 Hz, which leads to different weight and payloads. In the present study, the weight of FMAV and the extra moment supplied by horizontal tail are all tunable as the occasion demands in order to satisfy the force and moment equilibrium.

The cruising speed in forward flight is obtained based on the force equilibrium (averaged thrust equals averaged drag). According to the calculated results using CFD and FSI method, the average lift forces of FMAVs with rigid wing and flexible wing are 0.392 N and 2.734 N, respectively. Thus, the weight of FMAVs is 40 g and 279 g corresponding different wing according to the equilibrium relationship that weight equals lift forces. And the angle of attack is not treated as a known input parameter for the FMAV with rigid wing until force balance is satisfied when the angle of attack is set to 4°, as shown in Figure 5. As is seen in Figures 5 and 6, conditions in equilibrium flight are all satisfied when the angle of attack equals 4° and velocity is 4 m  $\text{s}-1$ and 14 m  $\text{s}-1$ , respectively, for FMAV with rigid wing and flexible wing. Table 4 presents the flow properties of trimmed computational model.

Figure 5.

Determination of the equilibrium between average thrust and drag of FMAV with rigid wing by varying velocity and angle of attack.

Figure 6.

Determination of the equilibrium between average thrust and drag of FMAV with flexible wing by varying velocity.

Table 4.

Computational model of trimmed FMAV with different wing.

Trimmed FMAV	Flexible-wing	Rigid wing
Angle of attack	4°	4°
Flapping amplitude	70°	70°
Flow velocity	14 m/s	4 m/s
Air density	1.2 kg/m3	1.2 kg/m3
Reynolds number	8.42 $\times 104$	2.41 $\times 104$
Flapping frequency	10 Hz	10 Hz
Reduced frequency	0.20	0.69
mass	0.392 N	2.734 N
Pitching moment of inertia $I_{\mathrm{yy}}$	0.0009163  kg ċ m2	0.0009163  kg ċ m2
Reference chord	0.088 m	0.088 m

Three complete flapping cycles have been calculated and results for the third cycle are only displayed. Figures 7 and 8 show the force coefficient of FMAV with rigid wing and flexible wing, respectively, when they are in the equilibrium state.

Figure 7.

Force coefficient of FMAV with rigid wing during third flapping period.

Figure 8.

Force coefficient of FMAV with flexible wing during third flapping period.

Figure 9 presents the snapshots of deformation during a half flapping cycle when the wing is flapping down and the angle of attack is 4°, the flow velocity is 14 m/s, flapping frequency is 10 Hz, and the flapping amplitude is 70°. The solid line profile represents the rigid flapping position, and the color contour surface corresponds the deformable wing. Figure 10 displays the surface pressure (relative to the free-stream pressure) at the mid-plane during down stroke for the flexible wing.

Figure 9.

Structural deformation of flexible wing during down stroke extracted from the CSD results.

Figure 10.

Distribution of surface pressure coefficients at the mid-plane during down stroke of flexible wing (top view on left, bottom view on right).

Results and analysis

Results of aerodynamic derivatives

After the equilibrium flight conditions have been determined, aerodynamic forces and moments are computed, respectively, when u and α vary independently and slightly around the equilibrium value. After subtracted equilibrium value from each quantity, the variation of the non-dimensional aerodynamic forces and moments against each of u and α is shown in Figures 11 to 14. As can be seen from the figures, X, Z, and M vary linearly against u and α, while $0\mathrm{deg}\le \alpha \le 8\mathrm{deg}$ and $0.9u_e\le \text{u}\le 1.1u_e$ indicate that there is a high linearization of the equations of motion and small disturbance assumption is reasonable. The non-dimensional aerodynamic derivatives X_u, Y_u, Z_u, $X_{\alpha }$ , $Y_{\alpha }$ , and $Z_{\alpha }$ are estimated using the data in Figures 11 to 14.

Figure 11.

The α-series force and moment data of FMAV with rigid wing.

Figure 12.

The u-series force and moment data of FMAV with rigid wing.

Figure 13.

The α-series force and moment data of FMAV with flexible wing.

Figure 14.

The u-series force and moment data of FMAV with flexible wing.

Since we only have $X_{\alpha }{Y}_{\alpha }{Z}_{\alpha }$ , here we take the same method as Taylor and Thomson used ³ to obtain $X_w{Y}_w{Z}_w$ . Hence, the w-derivatives are simply

$X_w=\frac{\partial X}{\partial w}\approx u_e^{-1}\frac{\partial X}{\partial \alpha }{Z}_w=\frac{\partial Z}{\partial w}\approx u_e^{-1}\frac{\partial Z}{\partial \alpha }{M}_w=\frac{\partial M}{\partial w}\approx u_e^{-1}\frac{\partial M}{\partial \alpha }$

(13)

Based on equation (13), all six static stability derivatives for different FMAV can be obtained and the system matrix in equation (2) can also be fulfilled. Here, we transfer the non-dimensional coefficients to dimensional form for convenience, which are shown in Table 5. Comparing aerodynamic derivatives between rigid wing and flexible wing, there is a large scale of difference, which means the flexibility of wing affects the aerodynamics derivatives and properties significantly. It seems to be that the flexibility makes the aerodynamic derivatives smaller in general and less sensitive especially for the Z_u, but except for the increase of the M_w.

Table 5.

Dimensional aerodynamic derivatives with respect to forward speed U and attack angle α for two FMAVs.

FMAV	$X_{\text{u}}$ $\mathrm{Nm}-1S$	$Z_{\text{u}}$ $\mathrm{Nm}-1S$	$M_{\text{u}}$ $N·S$	$X_{\text{w}}$ $\mathrm{Nm}-1S$	$Z_{\text{w}}$ $\mathrm{Nm}-1S$	$M_{\text{w}}$ $N·S$
Rigid	−0.564	1.986	1.606	−2.618	−20.759	−29.43
Flexible	−0.256	−0.082	0.544	−0.378	−4.733	−35.72

Solution of the small disturbance equations

The components in the system matrix of equation (2) are all acquired after the aerodynamic derivatives are calculated. The dynamic flight stability of the FMAVs can be analyzed by solving equation (2). The method for such a system of linear equations is presented in detail in many textbooks of flight dynamics. The results of eigenvalues and the corresponding eigenvectors are calculated and given in Tables 6 to 11.

Table 6.

Eigenvalues of the system matrix of FMAV with rigid wing.

Rigid wing	Model 1 ${\lambda }_{1,2}$	Model 2 ${\lambda }_3$	Model 3 ${\lambda }_4$
Eigenvalue	−10.015 ± 3.266i	−2.269	0.977

Rigid wing

As seen in Table 6, for the FMAV with rigid wing, there are a pair of complex eigenvalues with a negative real part representing a stable oscillatory motion (mode 1), and one negative real eigenvalue representing a stable subsidence motion (mode 2) and one positive real eigenvalue representing an unstable divergence mode (mode 3). This is similar to Taylor and Thomson’s results. ³

As is known to us, the disturbed motion is a linear combination of natural modes. When disturbed from its balanced flight, the disturbed motion of FMAV with rigid wing is a linear superposition of an unstable divergence mode, a pair of stable oscillatory modes, and a stable subsidence mode. Because the growth of the disturbed motion is determined by the unstable divergence mode, the system is unstable.

The pair of complex conjugate roots $\lambda =n\pm i\omega$ represents a damped oscillatory mode which has relatively short period (T = 1.924 s). The damping ratio ξ represents the damping of the motion

$\xi =-n/(\omega 2+n2)\frac{1}{2}$

(14)

which for the FMAV with rigid wing is ξ = 0.95. Because when ξ equals 1, it represents critical damping that the transition from sinusoidal to exponential motion, which denotes that the oscillatory mode is heavily damped.

What’s more, the static damping of the motion can be quantified by the time taken for an oscillatory mode to halve in magnitude, where

$t_{\mathrm{half}}=\frac{\ln 2}{|n|}$

(15)

This yields times to halve that are approximately the same as a single period of oscillation. ( $t_{\mathrm{half}}$  = 0.069 s). Since the flapping frequency is 10 Hz (the wing-flapping period is T = 0.1 s), thus the period of the oscillation for FMAV with rigid wing is about 19 times of the wing-flapping period, and the starting value of the oscillation will be half in less than one wing-beats.

$t_{\text{half}}$ or $t_{\text{double}}$ for the three modes and the period of the oscillatory mode and are presented in Table 7.

Table 7.

Dimensional time constants of the natural modes of FMAV with rigid wing.

Model 1			Model 2		Model 3
Stability	T	$t_{\text{half}}$	Stability	$t_{\text{half}}$	Stability	$t_{\text{double}}$
Stable	1.924 s	0.069 s	Stable	0.305 s	Stable	0.709 s

It is known to us that the magnitudes and phases of the disturbance quantities of each mode are determined by the corresponding eigenvector. In order to compare eigenvector quantitatively, we have scaled them to make $\delta \theta$  = 1 and in polar form (as shown in Table 8). Scaling the eigenvectors makes the comparison of the two different FMAVs more conveniently. $\delta w$ can be converted to radian unit by normalizing $\delta w$ by u_e in Table 8. Normalizing $\delta u$ by u_e similarly brings $\delta u$ into the dimensionless form as $\delta w$ .

Table 8.

Magnitudes and phase angles of the components of each of the three eigenvectors of FMAV with rigid wing.

Stiff	Model 1	Model 2	Model 3
$\delta u$ /uE	0.16(91.9°)	1.50(0°)	1.44(180°)
$\delta w$ /uE	0.95(143.5°)	0.04(0°)	0.09(180°)
$\delta q$	10.53(162°)	2.27(180°)	0.98(0°)
$\delta \theta$	1.00	1.00	1.000

We can see from Table 8 that for the oscillatory mode (first mode), the changes of forward velocity component ( $\delta u$ ) are insignificant. On the contrary, oscillatory changes in pitch rate ( $\delta q$ ) are notable. It seems to identify this oscillatory mode of FMAV with the short period pitching mode, which is intrinsically a rapid pitch oscillation with small varies in forward velocity.^28,29

On the contrary, change in the forward velocity component ( $\delta u$ ) and pitch rate ( $\delta q$ ) is obvious in the subsidence and divergence mode. But the phase of pitch rate ( $\delta q$ ) is completely opposite in these modes.

Flexible wing

Similarly, the eigenvalues and eigenvector of FMAV with flexible wing can be obtained using same method, as shown in Tables 9 to 11.

Table 9.

Eigenvalues of the system matrix.

Flexible	Model 1 ${\lambda }_{1,2}$	Model 2 ${\lambda }_{3,4}$
Eigenvalue	−2.359 ± 22.232i	−0.136 ± 0.300i

As can be seen from Table 9, there are two pairs of complex eigenvalues with a negative real part, representing two stable oscillatory motion. Different from FMAV with rigid wing, the flight dynamic stability of FMAV with flexible wing is inherently stable. And it is shown in Table 10 that there is a short period mode (T = 0.283 S) and long period stable mode (T = 20.94 s). It is noticeable that such flight stability is very similar to that of fixed-wing aircraft whose longitudinal equations of motion also have two pairs of complex conjugate roots, the short period motion, and a much longer period known as phugoid mode. The latter mode is an interchange between kinetic and potential energy, which resulted in a steady rise and fall in altitude. ³⁰ Normally, it is accompanied by slight changes in pitch and the angle of attack remains constant during the process. Similarly, when there is a disturbance of $\delta \theta$ for the long period mode of FMAV with flexible wing, variations in the forward velocity component ( $\delta u$ ) are notable in the oscillatory mode, as shown in Table 11.

Table 10.

Dimensional time constants of the natural modes.

Model 1			Model 2
Stability	T	$t_{\text{half}}$	Stability	T	$t_{\text{half}}$
Unstable	0.283	0.294	Stable	20.94	5.097

Table 11.

Magnitudes and phase angles of the components of each of the three eigenvectors.

Flexible	Model 1	Model 2
$\delta u$ /ue	0.048(88.9°)	2.154(112.8°)
$\delta w$ /ue	1.0(12.2°)	0.033(112.4°)
$\delta q$	22.357(96.1°)	0.329(114.3°)
$\delta \theta$	1.000	1.000

Discussion

The influence of the flexibility of wing of FMAV

By comparing the flight stability of FMAV with different wing, it is concluded that the flexibility can indeed change the flight stability properties from unstable to stable. And in real forward flight, as the present research subject, the Dove FMAV with flexible wing performs well and inherently stably, which proves the above simulated results reliable and correct. However, this may not be suitable for all FMAV. Because the FMAV in the present paper is special, it cannot represent all kinds of FMAVs. It can only prove that at certain case, the flexibility indeed improves the flight stability even from unstable to stable. In order to figure out whether other factor influences the flight stability, especially in the course of design for FMAV, here, the influence of position of mass center is studied carefully.

The influence of the position of mass center of FMAV

In the present part, the influence of the position of mass center of FMAV on the dynamic stability of FMAV with flexible wing is studied and analyzed. Here, we calculate the eigenvalues corresponding to mass center position l of FMAV from 0.05 to 0.3 in 1000 intervals (the reference length from position of mass center to leading edge of the wing) as shown in Figure 15. Table 12 shows the aerodynamic derivatives of two limiting cases: $l=0.05$ and $l=0.3$ . It can be seen from the table that the position of mass center only influences the moment related derivatives and $M_{\text{w}}$ changes much more than $M_{\text{u}}$ , which means that $M_{\text{w}}$ is more sensitive than $M_{\text{u}}$ .

Figure 15.

Shift of center of gravity(c.g).

Table 12.

Dimensional aerodynamic derivatives.

l	$X_{\text{u}}$ $\mathrm{Nm}-1S$	$Z_{\text{u}}$ $\mathrm{Nm}-1S$	$M_{\text{u}}$ $N·S$	$X_{\text{w}}$ $\mathrm{Nm}-1S$	$Z_{\text{w}}$ $\mathrm{Nm}-1S$	$M_{\text{w}}$ $N·S$
0.05	−0.256	−0.082	0.544	−0.378	−4.733	−35.72
0.3	−0.256	−0.082	1.306	−0.378	−4.733	0.178

The eigenvalues of system matrix of two limiting cases are shown in Tables 13 and 14. We can see that the mode changes from stable to unstable. In order to figure out how exactly the transformation develops, we calculate the results of system responding different $M_{\text{u}}$ and $M_{\text{w}}$ varying from stable state ( $l=0.05$ ) to unstable state ( $l=0.3$ ) in 1000 intervals as aforementioned. Since the eigenvalue is almost same when the $M_{\text{u}}$ changes within l changing, here we do not present the corresponding results.

Table 13.

Eigenvalues of the system matrix.

l	Model 1 ${\lambda }_{1,2}$	Model 2 ${\lambda }_{3,4}$
0.05	−2.359± 22.232i	−0.136 ± 0.300i

Table 14.

Eigenvalues of the system matrix.

l	Model 1 ${\lambda }_{1,2}$	Model 2 ${\lambda }_3$	Model 3 ${\lambda }_4$
0.3	1.186±2.119i	−1.8654	−5.496

The evolution of the root locus plots of the eigenvalues is plotted in Figures 16 and 17. As we can see from the corresponding figures, the dynamic stability changes gently at first and keeps the real part of eigenvalue negative. With the position of mass center changing, the first pair of complex eigenvalue turns splits to be two positive eigenvalue and the second pair of complex eigenvalue turns to be two negative value, which means the states changes from unstable to stable. During the whole course of transformation, most of states are stable and there is a critical position of mass center that determine the boundary of stable and unstable. If we can find the critical position of mass center, it will give an inspiration to the design of FMAV.

Figure 16.

Root locus plots of first pair complex eigenvalue.

Figure 17.

Root locus plots of second pair complex eigenvalue.

Conclusion

In the present paper, the effect of wing flexibility on the dynamic flight stability of FMAV in forward flight has been studied. An integrated aeroelasticity model of flapping wing flight was used to calculate the aerodynamic derivatives. The techniques of eigenvalue and eigenvector analysis were used to solve the equations of motion and analyze the stability properties. To study the influence of wing flexibility on the stability of flapping flight, flight stability of two different FMAV with flexible and rigid wing was modeled and compared. For the rigid wing FMAV, the system is unstable. But when the wing is flexible, it was found that the inherent stability of the system changed from unstable to stable. Even though this example may be an exception, it proved to be an example that wing flexibility can change or even improve the stability of FMAV. What’s more, the position of mass center is also a key factor for the stability properties of FMAV and could be applied in the design of FMAV.