Aerodynamics and dynamic stability of micro-air-vehicle with four flapping wings in hovering flight

Recently, a novel concept of flapping Micro-Air-Vehicles (FMAVs) with four wings has been proposed, which potentially utilizes the clap-and-fling effect for lift enhancement and agile maneuvers through an adjustment of wing kinematics. However, the application of the clap-and-fling effect in the four-winged FMAVs is underexplored and the dynamic stability is still unclear. In this paper, aerodynamics and flight dynamic stability of the four-winged FMAVs are studied experimentally and numerically. Results show that the clap-and-fling effect is observed when the flapping frequency is above 18 Hz. Due to the clap-and-fling effect, the lift generation and aerodynamic efficiency are both improved, which is mainly attributed to the fling phase. Further studies show that the clap-and-fling effect becomes weaker as the wing root spacing increases and is almost absent at a wing root spacing of 1.73 chord length. In addition, a wing with an aspect ratio of 3 can increase both lift generation and efficiency due to the clap-and-fling effect. Finally, according to the dynamic stability analysis of the four-winged FMAV, the divergence speed of the lateral oscillation mode is about 4 times faster than that of the longitudinal oscillation mode. Our results can provide guidance on the design and control of four-winged FMAVs.

Now, a series of high-lift mechanisms including delayed stall [6], wake capture [7], and clap-and-flip effect [8][9][10][11][12][13] have been clearly addressed. Flying stability and the controlling strategies used by insects have also been systematically revealed. With the theories how insects fly in mind, scientists turned their focus into how to fully apply those flying principles to a man-made FMAV and conducted systematical researches on aerodynamic, structure, power and control problems. For example, extensive experimental and numerical studies have been conducted on the effects of wing geometrical [14][15][16] (including wing aspect ratio, area, span et al.) and on both lift and power consumption. Wing passive deformations beneficial for lift augment and power efficiency enhancement have also been discovered [17][18][19][20]. Then, various MAV prototypes using gear transmission, piezoelectric as well as electromagnetic drives were built, such as Robobee [21] of Harvard University, hummingbird [22] of AeroVironment Inc., KU Beetle [23] of Konkuk University, and Colibri [24] of ULB. After flying test of both vertical take-off and hovering flight, the dynamic stability of those bionic FMAVs has been further studied to guide flight control system design so as to finally achieve maneuvering flight.
Although various FMAV prototypes are capable of free flight, rigorous requirements in size and weight not only impose great difficulty on MAV conceptual design and subsystem design but also limit MAV's performance. Lift enhancement is always desired to carry more payloads like camera to conduct specific tasks. In addition, traditional FMAV consists of two wings which are not only used for lift and thrust generation but also taken as controlling surfaces to produce torques for maneuvering flight, which brings great difficulty in accurately controlling the flight of FMAV.
In terms to enhance lift production and reduce difficulty in the design of MAV control system, four-winged FMAV combining the concept of FMAV and rotary wing has been proposed, such as Delfly Nimble [25] and RoboBee X-wing [26] (Fig.1). A fourwinged FMAV consists of two pairs of flapping wings. Lift augment was designed from the increase in the number of wings, and also benefited from the potential clap-andfling effect of adjacent wings at the end of both upstroke and downstroke. The clapand-fling effect was firstly proposed by Weis-Fogh [8] to explain how small wasp generated enough lift. He pointed that wing clap-and-fling effect will induce a circulation around the wing which accounts for the high lift production. This assumption was later confirmed by many scientists [9][10][11][12][13]. Although this principle was generally found used by natural tiny insects whose Reynolds (Re) number is lower than 100, researchers still explored to introduce it into MAV design. Hoang Vu Phan designed a FMAV model Fig. 1 Typical four-winged FMAV. a Delfly Nimble, b RoboBee X-wing with a flapping amplitude of 180°and wings clapped and flinged at the end of upstroke and downstroke. They measured the averaged lift of this model and found that the clap-and-fling effect leads to a lift augment by about 16% [23]. In summary, wing that could benefit from clap-and-fling effect in lift augment has been widely accepted, however, its effect on wing aerodynamic power consumption as well as power efficiency remains unclear and detailed investigation is still needed to design a micro four-wing FMAV for larger lift production.
With the increase in the number of flapping wings, it is possible for the four-winged FMAV to achieve maneuver flight by controlling the overall motion of the flapping wings like that used by multi-rotor drone. Even so, the dynamic stability characteristics of this type of FMAV should still be clarified firstly. In recent years, dynamic flight stability of insects has been revealed and it was found that insect hovering flight is dynamically unstable [27,28] to enhance maneuverability. Same with insects, bionic FMAV is also of the same dynamic stability property. The four-winged FMAV combines the concepts of FMAV and rotary wing. Unlike FMAVs, multi-rotor drone is capable to stabilize its flight. Thus, it raises the question that whether the flight of the micro four-wing FMAV is inherently stable, which must be answered before the design of FMAV controlling system.
In this paper, aerodynamics and flight dynamic stability of the four-winged FMAV was studied experimentally and numerically. The paper is organized as follows: firstly, both experimental as well as numerical methods (including computational fluid dynamic (CFD) simulation and dynamic stability analysis) are introduced. Then, based on a typical model, to what extent the clap-and-fling effect could affect lift, power consumption as well as power efficiency are studied. Effect of wing geometrical parameters, i.e. wing root distance and aspect ratio, on aerodynamic performance has been further studied using CFD methods. Finally, the longitudinal and lateral dynamic stability of the four-winged FMAV model are analyzed.

Lift-measurement platform and experimental model
Averaged lift (L) of the four-winged FMAV model was measured experimentally using a fast lift measurement method to reveal the clap-and-fling effect on averaged lift production, as shown in Fig. 2. The flapping mechanism was fixedly connected to the optical plate through the support column. The regulated power supply provided steady electricity for the flapping mechanism. The laser frequency measurement (measurement accuracy is 0.01HZ) was fixedly connected above the wing and the laser shed onto the top surface of the flapping mechanism at wing root. When the wing flapped and went through the laser, the laser was reflected back by the white resin, and the digital counter was increased. By counting how many times the wing flapped under the measurement during one second, we computed the flapping frequency. A precision electronic scale with measurement accuracy of 0.01 g was used for averaged lift measurement. It was not directly connected with the FMAV model and was just placed below the wing, and the distance between the wing and the electronic scale was about one span. When the wing flapped, the downwash air flowed to the plate. Therefore, the force measured by the precision electronic scale was equal to the lift generated by the flapping wing. It has been tested that this method measures lift with an error of about 10% [22]. Two pairs of flapping wings of the four-winged FMAV were completely symmetrical in structure and function, so only one pair of flapping wings was investigated in the experiment. The 3D printed FMAV model (Fig. 2b) was driven by a 7.4 V hollow cup motor, which was decelerated by a gear set with the reduction ratio of 20. Wing flapping amplitude (ϕ f ) of the model was about 90°and the maximum flapping frequency (n) was about 30 Hz. The wing (Fig. 2c) was made of a leading-edge carbon rod (1 mm diameter), a reinforcing carbon rod (0.5 mm diameter carbon rod) and a vertical carbon rod (1 mm diameter carbon rod) covered by a polyimide film with the thickness of 0.025 mm. The span of wing was 9 cm and the chord length of the wing was 4 cm. During wing flapping, the wing membrane was deformed with the acting of rod inertia force and aerodynamic pressure to achieve angle of attack. Using a high-speed camera, the experiment captured the motion of the wing at the end of the upstroke (Fig. 2d) when the flapping frequency was 25 Hz, at which time a significant clap-and-fling motion can be observed. In addition to measuring the lift of the paired wings, we also measured the lift of a single wing ( L S ), which was not contributed to clap-and-fling effect. Comparing 2 L S and L at the same flapping frequency, the clap-and-fling effect on averaged lift can be revealed.

CFD method
The experimental measurement can only obtain the averaged lift and it is difficult to analyze the clap-and-fling effect on lift, energy consumption and aerodynamic efficiency. Therefore, a simplified model was built according to the experimental test model and  using CFD method to comprehensively explore the influence of clap-and-fling motion on the aerodynamic performance (lift, energy consumption and lift efficiency) of the four-winged FMAV.
The wing was simplified to a rectangular wing that flaps horizontally. The geometry of the wing was determined by aspect ration AR (defined as the ratio of the length R to the chord length c), the distance between the body center and the wing root C R (C R = 0.7c), and the root spacing (d) between the adjacent wing roots at the end of upstroke. In order to define wing kinematics, the inertial coordinate system o g -x g y g z g , the body-fixed coordinate system o b -x b y b z b , and the wing-fixed coordinate system o-xyz were introduced (Seen in Fig. 3). The o g -x g y g plane of the inertial coordinate system was horizontal. The origin of the body-fixed coordinate systems o b was located at the center of the body. The o b-x b y b plane of the body-fixed coordinate system coincided with wing flapping plane, and o b x b was symmetric axis of the body in vertical plane and pointed backward. The wing-fixed coordinate system o-xyz moved with wing flaps with oz axis normal to the flapping stroke plane and always pointed upward. ox was tangent to the flapping trajectory of the wing, pointing in the direction of downstroke. The flapping angle (ϕ) and angle of attack (α) were defined using a sinusoidal function as Eqs. 1 and 2, respectively.
where T was the non-dimensional flapping period, α was the angle of attack in the middle of upstroke, and Δτ r (Δτ r = 0.25 T) was time for wing performing a flipping motion. In the CFD simulation, ϕ f = 90°and α =40°. Refer to wing kinematics of those insects using clap-and-fling effect for high lift production, the pitch axis varied during the clap-and-fling motion. When the FMAV clapped its wings, the pitch axis was the leading edge and it turned into the trailing edge after wing entered fling phase. To eliminate the changes in wing vertical displacement caused by the changes in the pitch axis, an additional vertical velocity of wing was introduced as presented in Eq. 3.
Based on wing kinematics defined, Re was defined based on the velocity (U = 2ϕ f nr 2 ) at the second moment of wing area (r 2 ) of the flapping wing area as Re = cU/ν, where ν was viscosity coefficient. In this paper, Re was set at 3725 based on the experimental measurement of the four-winged FMAV model.
The commercial software FLUENT was used to calculate the aerodynamic forces and power consumption. Considering that the four flapping wings were symmetrically arranged along the longitudinal and lateral directions of the FMAV respectively, only one flapping wing was calculated and the influence of the adjacent two flapping wings was simulated by setting symmetric boundary conditions. The simulation used dynamic mesh technology to realize the unsteady motion of the wing. The wing surface was constructed by 80,000 structural grids with unstructured grid wrapped. To improve the grid quality, the unstructured grid zone was further divided into two parts: the inner part was dense unstructured grid with a radius of 20c around the wing surface and the outer was a sparse unstructured gird with a radius of 40c. (seen in Fig. 4).
The 3D pressure-based coupled solver with laminar model was used for the transient simulation of the flapping wing. We defined the lift normal to the wing flapping plane C D and C P were further calculated as C L ¼ L 0:5ρU 2 S ; C D ¼ D 0:5ρU 2 S ; C P ¼ P 0:5ρU 3 S , where S was wing area. We used C L , C P to denote the average lift coefficient and aerodynamic power coefficient during one stroke cycle, and C D was defined as average drag coefficient during half stroke cycle. Based on C L , C P , aerodynamic efficiency η was computed as η In order to test the independence of the grid used in the numerical simulation, we have verified the grid density, the boundary domain size and the time step separately. During grid-independence test, the parameters of the flapping wing were selected as Re = 3725, d = 0.28c and AR = 3. Under the conditions of 20 times far field and 800 steps/cycle, mesh 1 (220,000 grids), mesh 2 (360,000 grids), and mesh 3 (650,000 grids) were selected for grid density verification. The results are shown in Table 1. It can be seen from the comparison of the data that the calculation results of the medium density grid (mesh 2) and the dense grid (mesh 3) were relatively close, and the results of the sparse grid (mesh 1) calculation were quite different. Therefore, in order to ensure the calculation accuracy and improve the calculation efficiency at the same time, a 360,000 grid was selected for subsequent calculation. The same method was used to verify the boundary domain size and calculation time step. Finally, the numerical simulation of this study was carried out using a grid with a far field size of 20c and a time step of 0.005 T.
In addition, in order to verify the difference between the single-wing simulation and the full simulation, the study also compared the results of symmetric processing with that of full model as shown in Table 2. It can be seen from the table that the results obtained by the two methods were basically the same. Therefore, it was feasible to use the symmetric processing to perform the aerodynamic simulation of the multi-wing wing model.
To validate the CFD method used in this paper, CFD simulation toward a flapping wing with Re = 16,100, AR = 3, ϕ = 120°and α = 45°was conducted and compared with the results calculated by Hao et al. [30]. Instantaneous C L and C D in a flapping cycle from two methods were shown in Fig. 5. As can be seen from Fig. 5, instantaneous C L and C D from these two methods were approximately the same despite some slight differences in the reversal phase. Besides, the stroke-averaged aerodynamic coefficient calculated using the method in this paper was also nearly the same as that calculated by Hao et al. with an error less than 2%. In summary, the CFD simulation method used in this work was reasonable.

Dynamic stability model
Based on the idea of insect dynamic stability analysis [27,31], a six-degree-of-freedom rigid body dynamics equation was established based on the average aerodynamic force and moment (Fig. 3a). The longitudinal and lateral motion equations of the aircraft were as Eq. 4.
where, u, v, w, p, q, r were the corresponding projections of the centroid velocity and angular velocities of the FMAV on the x b , y b , z b axis, respectively. θ was the angle between the o b x b axis and the ox axis, and φ was the angle between the o b y b axis and the oy axis. I xx , I yy and I zz were the moment of inertia around the x b , y b , z b axis, respectively, whose parameters were from the experimental model (I xx = 5.634 × 10 − 5 kg·m 2 , I yy = 84.296 × 10 − 6 kg·m 2 , I zz = 4.810 × 10 − 5 kg·m 2 ). X, Y, Z, L, M, N were the corresponding projections of aerodynamic forces and moments on the x b , y b , z b axis. The small perturbation equation (shown as Eq.5) can be obtained by linearizing Eq. 4 in the equilibrium state and ignoring the small amount above the second order. where: Where m was the mass of the two pairs of FMAV and m = 62 g. Because the longitudinal motion equations and lateral motion equations can be decoupled, we dealed with these equations separately.
For the aerodynamic derivatives, we can get the accurate results of this parameters but at a high cost of calculation sources and calculation time, which did not satisfy requirements in the analysis of dynamic stability especially in conceptual design phase. Therefore, the quasi-steady aerodynamic estimation method was adopted to calculate the aerodynamic derivatives. A detailed description and verification of this method can be found in the literature of these researches [27,28,31].

Aerodynamic performance of a typical micro four-wing FMAV
Averaged lift of a typical micro four-winged FMAV model was firstly measured to verify the effect of the clap-and-fling effect on lift production. To get rid of the error coming from wing structure and experimental measurement, five pairs of wing are measured and the experiments are repeated three times for each trial. Averaged lift (L) produced by two wings with and without the clap-and-fling effect at different flapping frequencies is presented in Fig. 6. It can be seen from Fig. 6 that whether the clap-andfling effect could bring lift augment depends on wing flapping frequency. When the flapping frequency is lower than a critical value, which is 18 Hz in this experiment, there is no significant change in lift production. Once the flapping frequency is higher than 18 Hz, lift produced by wing with clap-and-fling effect motion becomes larger than that without the motion and the lift difference caused by clap-and-fling effect almost increases linearly with wing flapping frequency. In the studied range of flapping frequency, the lift is increased by 11.2% due to clap-and-fling effect when f = 25 Hz. As shown in Fig.2d, the clap-and-fling effect of two adjacent wings is mainly caused by wing's flexible deformation. At a flapping frequency smaller than 18 Hz, the flexible deformation of the flapping wing is not significant and as a result, clap-and-fling effect is also weak. As wing flapping frequency increases, wing passive deformation increases and there is obvious clap-and-fling motion of wing trailing edge at the end of upstroke, which consequently leads to a higher lift production at this phase.
To further test the clap-and-fling effect on lift, power consumption and aerodynamic efficiency at high flapping frequency, a CFD simulation is conducted. Note that for the experiment model, wing always pitches around its leading edge. Differently, for the model used in CFD simulation, wing pitching axis during the reverse phase is changed between leading edge and trailing edge, which resembles that of tiny insect based on the observation results, to fully utilize the clap-and-fling effect. Parameters used in this simulation are as follows: Re = 3725, d = 0.28c, AR = 3, ϕ = 90°. C L , C D , C P and η computed are presented in Table 3. Table 3, clap-and-fling effect leads to wing C L increasing by 19% from 2.044 to 2.436 at a cost of a larger power consumption, which increases by 23.9%. As the increase in lift is more pronounced, η also increases by about 4.8%. It suggests that wing clap-and-fling effect is not only beneficial for lift augment but for power efficiency enchantment. Figure 7a plots the instantaneous aerodynamic forces and the power consumption in a flapping period. As can be seen from Fig. 7, C L has a periodic change over time with a period at 0.5 T, while C D seems changes without period which is caused by the definition of L and D. In fact, C D in upstroke and downstroke has the same absolute value while the sign of C D is changed. For both wings with or without clap-and-fling effect, there is a larger contribution of aerodynamics observed due to the rapid pitching of the wing. The increase in lift, drag and power consumption caused by the clap-and-flip effect mainly occurs when wings fling (e.g. t = 0~0.125 T and 0.5~0.625 T). In the four-winged FMAV, the clap-andfling motion repeated twice for a specific wing at the end of both downstroke and upstroke. Consequently, there are two significant increment in instantaneous lift in a flapping period.

As shown in
In order to further analyze the clap-and-fling effect on the aerodynamic forces production at Re~ο (10 3 ), vorticity field around wing with and without clap-and-fling  Fig. 7d). In the fling phase (e.g. t = 0.06 T), a significantly stronger attachment leading edge vortex (LEV) and tip vortex (TV) is formed in the upper surface of the wing due to a rapid pitching. The zone occupied by the LEV and TV is of low pressure, which results in a large normal suction acting on the wing leading to large lift and drag. As seen in Fig. 7d, at this time, LEV and TV formed around the wing with clap-and-fling motion is much intensive with larger strength, so there is much larger lift and drag produced and instantaneous C L , C D and C P increases by 48.3%, 53.6% and 14.5%, respectively. When the wing translates (e.g. t = 0.32 T), the fling motion of wings has already finished and the clapand-fling effect becomes so weak and the instantaneous structure of the vorticity field around the wing with or without the clap-and-fling effect is basically the same, so the instantaneous aerodynamic force is basically the same at this time. When the wings move to the end of upstroke and the clap motion begins (e.g. t = 0.42 T), the LEV remains attached to the lower surface of the wing. Although a downwash flow is induced, the intensity of the LEV is not influenced significantly. As a result, the wing aerodynamics and power consumption of wing with clap-and-fling effect motion slightly Table 3 Average aerodynamic coefficient and efficiency of wing with or without clap-and-fling effect  changed and C L , C D as well as C P increase by 12%, 15% and 6%, respectively, compared to its counterpart.

Influence of geometric parameters on aerodynamics
This section explored the effect of wing geometrical parameters on the aerodynamic performance enhancement led by clap-and-fling effect. Influence of wing root spacing and aspect ratio of the wing was studied using CFD simulation.

Wing root spacing
Three wing roots spacing, that is d = 0.28c, 0.43c, 1.73c, are considered with Re, AR, ϕ unchanged. C L , C D , C P and η computed at different d are shown in Table 4. C L , C D and C P are plotted in Fig.8. Results of wing without clap-and-cling effect are also presented for comparison. It can be seen from Table 4 that within the studied range of d, wing lift and power efficiency could always be improved by the clap-and-fling effect.
When the d is 0.28c, increase in C L can be maximized and its value increases from 2.044 to 2.503 (by 23%). Despite a significant increase in C P (increased by 19%), aerodynamic efficiency in lift production is still slightly increased by 7.4%. As wing root spacing increases, clap-and-fling motion between two adjacent wings is weakened, and lift augment and power efficiency caused by the clap-and-fling effect are all gradually reduced. When the wing root spacing increases larger than 1.73c, the aerodynamic performance of these two types of flapping wing is basically the same, which means the clap-and-fling effect is negligible. As can be seen from Fig. 8, C L , C D and C P produced during two reversal phases reduced with the increment of d (seen from t = 0.06 T, 0.42 T, 0.56 T and 0.92 T), especially in the fling phase (e.g. t = 0.06 T). At t = 0.06 T, wing lift coefficient of d = 0.28c, 0.43c and 0.73c is increased by 55.5%, 32.7% and 3.3% compared to wing without clap-and-fling effect, while it only increased by 13.2%, 8.3% and 0.5% at t = 0.42 T (i.e. clap phase). In summary, a smaller wing root spacing is suggested to fully utilize the "clap-and-fling" effect for a better aerodynamic performance.

Aspect ratio
As the aspect ratio of insect wings is in the range of 1~5, therefore three typical aspect ratio, that is AR = 1, 3 and 5, are considered. Please note that the chord length is unchanged with aspect ratio. To remain Re = 3725, the non-dimensional flapping frequency is changed with aspect ratio. Wing with small aspect ratio flaps fast than that with large one. C L , C D , C P and η of wing with different aspect ratio are shown in Table 5 while time courses of C L , C D and C P in a flapping period are presented in Fig. 9. As can be seen from Table 5 that at a small AR (about 1), the clap-and-fling effect plays a reverse role in lift production while power consumption is slightly enlarged. C L of wing with clap-and-fling effect reduced from 6.145 to 5.952 by 3.4% compared to single wing. At a constant Re, wing with small aspect ratio flaps at high frequency and larger aerodynamic contribution is observed in fast reversal phase at the end of both upstroke and downstroke. Compared to single wing, increment in instantaneous lift of wing with clap-and-fling effect is greatly reduced compared to cases with AR > 1. After wings fling and undergo translation, loss in lift further offsets the lift augment from clap-and-fling effect. Taking a complete flapping period into account, wing's averaged lift and power efficiency decreases. As AR becomes larger than 1, lift augment and power efficiency can be got by clapand-fling motion. Furthermore, compared to case at AR = 5, wing with an aspect ratio of 3 can benefit more in lift augment. As seen from Fig. 9, the instantaneous aerodynamic forces of wings with or without clap-and-fling effect at reversal phase both decreased significantly as wing flapping velocity decreased with large AR. From instantaneous aerodynamic coefficient and power coefficient shown in Fig. 8, the peak value of C L is found about 1.5 times larger than that produced by the single wing when AR ≥ 3, causing lift augment in those cases at AR > 1. Further comparison of the instantaneous aerodynamic forces produced by wing with AR = 3 and AR = 5 shows that in case of a large aspect ratio, there is a small amount of lift loss (about 5.8%) occurred during the translation phase compared to single wing while there is a lift increment (about 6%) occurred during the translation when AR = 3, so the overall lift augment of large aspect ratio wing is slightly smaller than that of wing with medium aspect ratio. To summarize, a medium aspect ratio wing is more recommended to design a micro four-wing FMAV when Re is kept constant.

Dynamic stability of a four-winged FMAV under typical design conditions
The dynamic stability analysis is based on the small perturbations from equilibrium state, so it is necessary to estimate the kinematic parameters of the FMAV in equilibrium state. Here, the quasi-steady aerodynamic method is applied to achieve the equilibrium state of aerodynamic forces and moment. The equilibrium state parameters are roughly as follows: α =45°, ϕ f = 90°. In order to solve the small perturbation equation, the stability derivatives need to be calculated firstly in both longitudinal and lateral equations. The calculation method is the same as that used by Sun and Xiong [27] and Sun et al. [28]. Here is an example of the longitudinal U-series aerodynamic derivatives (X u , Z u and M u ): it is necessary to take a number of effective small disturbances (usually within 10% of the reference velocity), and set w and q to zero, then calculate the aerodynamic forces and moments curve with u. The slope of the curve at the equilibrium state is thus the corresponding aerodynamic derivative. The other aerodynamic derivatives in longitudinal and lateral are computed similarly. Finally, the non-dimensional aerodynamic derivatives of the four-winged FMAV can be obtained as shown in Table 6.
Bring the non-dimensional aerodynamic derivative into the system matrix, we can get longitudinal and lateral equations as Eq. 6 and 7, respectively.  Then the techniques of eigenvalue and eigenvector analysis are applied to solve the longitudinal equations and the lateral equations of the four-winged FMAV, and the corresponding eigenvalues of longitudinal stability matrix A 0 and lateral stability matrix A 1 are shown in Table 7, respectively. It can be seen that the longitudinal mode and the lateral mode have three modes respectively: mode 1 (unstable oscillation mode), mode 2 (slow subsidence mode) and mode 3 (fast subsidence mode).
The period of longitudinal mode 1(T osc ) and the time to double the initial perturbation (t double ) (or half the initial perturbations, t half ) can be calculated by Eq. 8.
where,n andω are the real and imaginary parts of the eigenvalue, respectively. Thus, the time constants (dimensionless by the stroke period) of the above modes are shown in Table 8. For the oscillating divergent modes of the FMAV, the longitudinal perturbation needs 93.65 strokes to diverge twice, while the lateral perturbation needs only 21.3 strokes to diverge twice. So, the lateral control response of the four-winged FMAV needs to be faster. The eigenvector of each mode determines the magnitude and phase relationship of state variables. The eigenvectors of the longitudinal stability matrix A 0 and the lateral stability matrix A 1 are shown in Table 9. The longitudinal and lateral eigenvectors are normalized by δθ and δφ, respectively. As can be seen from Table 9, the longitudinal mode 1 is dominated by horizontal forward motion (u) and pitch motion (θ and q), and the vertical motion is about one order smaller in magnitude and can be ignored in this mode. Since the phase difference of u and q is less than 90°, the motion is horizontal forward motion with nose-up and horizontal backward motion accompanied by nosedown. Longitudinal mode 2 is dominated by horizontal forward motion and pitch motion, too. But the phase difference between u and q is 180°, indicating that there is a nose-down angular velocity in nose-up state and an nose-up angular velocity in the nose-down state, which means there will be a tendency to overcome the disturbance and the movement tends to be stable. Longitudinal mode 3 is dominated by ascending and sinking motions. The lateral mode 1 is dominated by the side slip motion (v) and the roll motion (p), and the yaw motion is of a smaller magnitude. Since the phase difference of v and p is greater than 90 degrees, the motion is a right horizontal movement with rolling left or left horizontal movement with rolling right. The lateral mode 2 is dominated by yaw motion. The lateral flight mode 3 is dominated by the side Table 6 Non-dimensional aerodynamic derivatives motion (v) and the roll motion (p), too. But the phases of δv+, δp + and δϕ + indicate that the motion is a translation motion to left with a right roll angle accompanied by rolling left. Rolling angular velocity is opposite to the roll angle.

Discussion
In order to further explore the dynamic stability influence of the four-winged FMAV, we compared both the longitudinal mode and lateral mode of four-winged FMAV with those of some insects and a pair of flapping wings MAV (as shown in Table 10 and Table 11, respectively) [28,29,31]. It can be found that they have the same longitudinal modes. The longitudinal modes of the four-winged FMAV, insects and a pair of flapping wings are divided into three types: the unstable oscillatory divergence mode, the stable fast subsidence mode and the stable slow subsidence mode. For the control design of aircraft, we are more concerned with the unstable oscillatory divergence mode. The time to double the initial perturbation is determined by the real part of eigenvalues. It can be inferred that the t double of four-winged FMAV is about one order larger in magnitude than that of the KU Beetle, thus, leaving enough margin for the response time of the control system. From the longitudinal eigenvalues, we can see the mode of four-winged FMAV is more like the mode of insects. Thus, the four-winged FMAV can have more possibility to maneuver as the natural insects. While the four-winged FMAV's lateral mode is not the same as some insects (i.e. BB, DF, HM). Insects like BB, DF, HM have one stable oscillatory mode, one stable subsidence mode and one unstable divergence mode. But the four-winged FMAV has one unstable oscillatory divergence mode and two stable subsidence mode, the same as its longitudinal mode. The modal structure of the four-winged FMAV is more appropriate to the modal structure of HB. Although the oscillatory subsidence mode of fourwinged FMAV is unstable, its time to double the initial perturbation is smaller than BB, DF and HM. Thus, we can call the unstable oscillatory divergence mode as a nearly neutrally stable oscillatory mode [29]. It means that we can control this unstable mode in time through the control system.
Considering the longitudinal and lateral mode comprehensively, the four-winged FMAV is more advantageous for the control system, because it provided enough response time for the control system.

Conclusion
A four-winged FMAV, which utilizes the clap-and-fling effect to improve its aerodynamic performance, is proposed in this research. Experimental measurements of wing kinematics are conducted and numerical simulations are employed to unfold the effect of clap-and-fling effect on the four-winged FMAV. A dynamic stability analysis is then carried out using the modal analysis method. Our findings are concluded as follows: 1. For a flapping wing with an aspect ratio of 3, the clap-and-fling effect is observed when the flapping frequency is above 18 Hz. Within the parameter space of our study, the clap-and-fling effect results a linear lift enhancement over the flapping frequency. Due to the clap-and-fling effect, the lift enhancement is around 11.2% at a flapping frequency of 25 Hz. 2. The effect of clap-and-fling effect on the lift enhancement of the four-winged FMAV is intensified as the wing root spacing decreases. Moreover, the lift enhancement is mainly attributed to the fling phase. The clap-and-fling effect is almost absent when the wing root spacing is above 1.73 chord length. 3. A further study on the aspect ratio (AR) shows that the lift enhancement due to the clap-and-fling effect is observed at AR = 3 or 5. In addition, the improvements of cycle-averaged lift generation and efficiency for a flapping wing with AR = 3 are 19.2% and 4.8%, respectively. However, as a result of the clap-and-fling effect, the cycle-averaged lift generation and efficiency of a flapping wing with AR = 1 is reduced by 3.1% and 5.6%, respectively. 4. The longitudinal and lateral dynamic stability modes of the four-winged FMAV are akin to insects and can be separated into three components: the unstable oscillatory divergence mode, the stable fast subsidence mode and the stable slow  Our study can thus contribute to the design of four-winged FMAVs to better exploit the clap-and-fling effect, thereby extending its flight endurance under a heavy payload. Our dynamic stability analysis can shed light on the control of four-winged FMAVs.