Experimental and Numerical Investigation on Dragonfly Wing and Body Motion during Voluntary Take-off

We present a detailed analysis of the voluntary take-off procedure of a dragonfly. The motions of the body and wings are recorded using two high-speed cameras at Beihang University. The experimental results show that the dragonfly becomes airborne after approximately one wingbeat and then leaves the ground. During this process, the maximum vertical acceleration could reach 20 m/s2. Evidence also shows that acceleration is generated only by the aerodynamic force induced by the flapping of wings. The dragonfly voluntary take-off procedure is divided into four phases with distinctive features. The variation in phase difference between the forewing and hindwing and angle of attack in the down-stroke are calculated to explain the different features of the four phases. In terms of the key parameters of flapping, the phase difference increases from approximately 0 to 110 degrees; the angle of attack in down-stroke reaches the maximum at first and then decreases in the following take-off procedure. Due to experimental limitations, 2-D simulations are conducted using the immersed boundary method. The results indicate that the phase difference and the angle of attack are highly correlated with the unsteady fluid field around the dragonfly’s wings and body, which determines the generation of aerodynamic forces.

Slow take-off. Chen et al. 14 examined the voluntary take-offs of droneflies, Eristalis tenax, and showed that droneflies do not jump but only flap the wings to lift themselves into the air. As a result, the take-off process takes a longer time.
In sum, the combined use of flapping and leg jumps results in a fast take-off process. However, it is difficult to design and manufacture a system using this mechanism for MAVs at the current stage.
Among the many insects with outstanding flying abilities, dragonflies possess unique flight characteristics and strong maneuvering abilities 15 . Experimental studies on the gliding, hovering and uniform forward flight of dragonflies have been reported [16][17][18] . In gliding flight, the wings of dragonflies remain almost immobile, and the forewing and the hindwing are almost on the same plane. In hovering and uniform forward flight, the dragonfly flaps the wings asynchronously with a phase difference between the forewing and hindwing of approximately 90° to 150°, and the positional angle and angle of attack of the wing are hardly different between each process. Sun and Lan 19 calculated the lift and thrust coefficient of dragonflies with different phase differences between the forewing and hindwing using a three-dimensional numerical simulation method. They showed that synchronized flapping can generate the largest lift and that the 90° phase difference between the forewing and hindwing is helpful for thrust. However, the take-off behavior of dragonflies has not been well studied, which is likely the most important and most powerful process. In this paper, the process of voluntary take-off of a dragonfly (Pantala flavescens) is recorded using two high-speed cameras and a three-dimensional reconstruction technique to study the real motion of the body and wings. Based on these results, the kinematics of the body and wings are calculated to illustrate the link between wing flapping and body movement, which results in an understanding of the voluntary take-off procedure of dragonflies.
Based on the take-off mechanism, the wingbeat of dragonflies could inform a mode of rapid take-off solely relying on flapping for MAVs. To further study why the wingbeat of a dragonfly can induce such high aerodynamic lift and to determine the key parameters related to the generation of aerodynamic lift, two-dimensional numerical simulation based on the IBM (immersed boundary method) is carried out for various flapping phase differences (between the forewing and hindwing) and angles of attack in the down-stroke.

Experimental Results and Discussions
Thirty dragonflies were tested. The maximum and minimum masses were 359 mg and 331 mg, respectively, with a relative difference of approximately 8%. The maximum and minimum body lengths were 49.4 mm and 47.9 mm, respectively, with a relative difference of approximately 3%. Due to these small differences in mass and size, the effects of mass and size were not taken into account.
During the experiment, parameters related to the voluntary take-off procedure were only recorded for eight individuals. Some similar phenomena were observed: First, dragonflies tended to leave the ground after one wingbeat. Second, dragonflies started to flap their wings synchronously, and the phase difference increased from 0 degrees during the subsequent wingbeats. Of the eight individuals, only three could be further analyzed because all the recorded parameters could be recognized only for these three individuals. In some cases, the wings were hidden behind other wings or the dragonfly flew outside the focused plane of the cameras.
For these three cases (labelled as D1, D2 and D3), all the analyzed results are presented, which is shown on Table 1.
The kinematic parameters share very similar values: In terms of the time of take-off procedure, for D1, D2 and D3, it is 350 ms, 380 ms and 300 ms, respectively. In terms of displacement, the largest differences are 0.08 m for upward movement and 0.02 m for forward movement. For average frequency, the largest differences are 0.7 Hz for forewing and 0.3 Hz for hindwing. For average stroke amplitude, the largest difference are 1.2 degree for forewing and 0.7 degree for hindwing. For the range of the phase difference, the minimum values are all 0 while the maximum fluctuating in the range from 109.2 degree to 110.6 degree. In the end, for the maximum value of the angle of attack in down-stroke, the largest differences are 1 degree for forewing and 4 degree for hindwing. Therefore, only the experimental results from one individual (D1) will be further analyzed to determine the mechanisms of dragonfly voluntary take-off.
Voluntary take-off process. There are approximately 8 stroke cycles in the voluntary take-off process of D1, which lasts for approximately 0.4 s. The sampling frequency of the camera is 1000 Hz, so 400 frames are recorded. After individual analysis of these 400 frames, 16 typical moments are chosen and displayed to describe the voluntary take-off procedure of the dragonfly (Fig. 1).
In Fig. 1, the 16 typical moments are labeled a to p. The time of moment a (corresponding to the beginning of wing flapping) is set to be 0. For the entire process, four phases with distinctive features can be identified: Phase I: from 0 to 47 ms. At 0 ms (moment a) the dragonfly starts to raise the fore and hind wings synchronously. Later, at 21 ms (moment c), both the forewing and hindwing reach the top position for the first time and start to flap down synchronously, while the body starts to rise up. Moment b (11 ms) is one point in the first up-stroke, and moment e (31 ms) is one point in the first down-stroke.
Phase II: from 47 ms to 80 ms. At 47 ms (moment f) the first down-stroke is almost over, and the dragonfly leaves the ground. From 47 to 80 ms, the dragonfly flaps both the fore wing and hind wing synchronously for Phase III: from 80 to 275 ms. At 80 ms, the flapping motion of the hind wing starts to become faster than that of the fore wing. As a result, the phase difference increases between the motion of the fore wing and hind wing. During the entire phase, the phase difference increases gradually from approximately 0 degrees to 110 degrees (the variation in phase difference is plotted in Fig. 2). Meanwhile, the rate of variation in the phase difference also changes; from 80 ms to 150 ms (moments j to l), the phase difference increases much more quickly than from 150 ms to 275 ms (moments l to n).
Phase IV: from 275 to 350 ms. In this phase (moments n to p), the phase difference is fixed at approximately 110 degrees. Finally, the take-off process is over, and the dragonfly starts to fly forward.
After three-dimensional reconstruction based on all recorded frames, the displacement of the body centroid in a relative coordinate system, which is detailed in the Methodology Section, is obtained and plotted in Fig. 3. Note that the displacement along the X-axis is very close to zero throughout the entire process of take-off, so only the displacements along the Y-axis and Z-axis are plotted. The coordinates at 0 ms are set to (0, 0, 0). All the typical moments consistent with Fig. 2 are also labeled.
In Phase I (moments a to f), only displacement in the upward direction (along the Z-axis) can be detected. In the following Phases II to IV (moments f to p), movements in both the upward direction (along the Z-axis) and the forward direction (along the Y-axis) are captured. The displacements in both directions increase at a nearly constant rate. As a result, from 0 to 350 ms, the dragonfly flies approximately 0.25 m upward and moves approximately 0.12 m forward during the entire process of voluntary take-off.

Body kinematics
To further study the voluntary take-off process, the body kinematics were calculated according to the displacements. The first derivative of the displacement-time function yields the velocity of the centroid, and the second derivative yields the acceleration of the centroid. The first and second derivatives of the displacements are  First, the variations in both velocity and acceleration correspond to different features of motion in the four phases. In Phase I, the dragonfly stays on the ground, so there is neither displacement nor acceleration along the Y-axis. In Phase II, both velocity and acceleration along the Y-axis begin to increase when the dragonfly starts to take off. However, the velocity and acceleration remain small in this phase. In Phase III, the velocities along the Z-axis and Y-axis keep increasing, but the accelerations along both directions finally decrease. This distinguished change may be induced by the change in the phase difference between the fore wing and hind wing, which will be further discussed in the simulation section of this study. In Phase IV, velocities and accelerations remain constant, corresponding to the constant phase difference.
Second, among the 8 peaks of acceleration along the Z-axis detected during this process, the first and second peaks labeled by A and B in Fig. 4b have the same largest value (20 m/s 2 , twice the acceleration of gravity). The first peak A results from the first down-stroke corresponding to moment d in Fig. 1, and the second peak B is induced by the second down-stroke corresponding to moment i in Fig. 1. Note that the acceleration curve does not show a sudden rise in the first two wingbeats. Thus, it can be verified that leg leaping is absent during the voluntary take-off of the dragonfly.
Wing kinematics. Figure 5 shows the measured angles of wing flapping as functions of time for the fore wing and hind wing. During the voluntary take-off of the dragonfly, the positional angle φ fluctuates, and the variation in the phase difference can be clearly seen. Let Φ and φ denote the stroke amplitude and mean stroke angle, respectively; they are defined by where φ max and φ min are the maximum and minimum values of φ, respectively. Figure 5 also shows that for both wings, the angles of attack (α) at the bottom of the down-stroke are approximately 85 degrees in Phase I and Phase II, much greater than those in the following phases. Actually, in Phase III, the angle of attack decreases from 85 degrees. Finally, in Phase IV, the angle of attack remains approximately 60 degrees.
To further analyze the wing flapping motions in Fig. 5, the frequency of wingbeat (n; estimated based on the time period between two adjacent peaks of the angle fluctuation), the amplitude of stroke (Φ) and the mean stroke angle (φ ) are calculated and summarized in Fig. 6. For the first two wingbeats, the estimated frequency is 17.2 Hz for the fore wing and 19.2 Hz for the hind wing. After that, the frequency fluctuates around 30 Hz in each wingbeat cycle.
The amplitudes of the stroke of both the fore wing and hind wing reach maximum values of 90.2 degrees and 92.0 degrees in the first wingbeat, respectively. Later, the stroke amplitude is almost constant. For the fore wing, the amplitude of stroke goes down until the 6 th wingbeat cycle and then increases. For the hind wing, the amplitude decreases to the minimum at the 4 th wingbeat cycle and then increases again. The mean stroke angles of both the fore wing and hind wing change from negative to positive values.
Summary of this section. The voluntary take-off of the dragonfly can be divided into four phases according to the distinctive features of body motion, body kinetics and wing motion. From Phases I to IV, the dragonfly is airborne in one wingbeat and then starts to take-off. The maximum vertical acceleration is close to 20 m/s 2 . The aforementioned experimental results show that the acceleration is generated by only the aerodynamic force. The phase difference between the forewing and hindwing increases from approximately 0 degrees to 110 degrees. The angle of attack (α) at the bottom point of the down-stroke of both wings is significantly greater at first and then decreases. During the four phases of the take-off procedure, the body motion of the dragonfly varies tremendously, indicating a variation in the aerodynamic forces generated by the flapping wings. This variation in force is accompanied by the change in two kinematic parameters of the wingbeat, namely, the phase difference and the angle of attack, as mentioned above. It would be of interest to investigate their influence on the aerodynamic forces.

Results of CFD (Computational Fluid Dynamics) simulations and discussions
Based on the experimental results, the kinematics of the voluntary take-off process of a dragonfly are described. In this section, two-dimensional immersed boundary computational fluid dynamics is used to investigate the key parameters of wingbeat related to the aerodynamic lift in the different phases of the take-off process.
A few studies have been done to assess the unsteady large-lift mechanisms of flapping insect wings [20][21][22][23] . When the insect wing is flapping, a leading-edge vortex (LEV) is formed, but the vortex does not detach from the wing in an entire down-stroke or up-stroke. Thus, a large aerodynamic force produced by the vortex can be maintained. Therefore, the large aerodynamic force in insect flight is closely related to such a vortex. This mechanism is also called the delayed-stall mechanism because the stall is delayed permanently. It is also known that the aerodynamic lift is produced mostly during the down-stroke, while the drag is produced mostly during the up-stroke.
In the experiments mentioned above, we found that the phase difference and angle of attack in the down-stroke vary during the voluntary take-off of the dragonfly. The influences of these two parameters on the vortex structure and the aerodynamic forces are analyzed with the two-dimensional immersed boundary method.  The lift coefficient (C L ) and thrust coefficient (C T ) of wings during one wingbeat are calculated and integrated to derive the cycle averaged lift coefficient (C L, avg ) and thrust coefficient (C T, avg ), as shown in Fig. 7. Note that positive C L represents lift and positive C T represents thrust. Figure 7a shows that both lift and thrust coefficients in case 1 fluctuate more heavily than those in case 2. For example, the peak lift coefficient in case 1 reaches 58, while it is only approximately 40 in case 2. The trough lift coefficient is approximately −40 in case 1 but only −20 in case 2. The average lift coefficient in case 1 reaches 4.5, 12.5% greater than that in case 2, but the average thrust coefficient in case 1 is smaller than that in case 2, as presented in Fig. 7b.
The values of the dimensionless pressure P at τ = 0.43 in Case 1 and Case 2 are shown in Fig. 8. The pressure differences between the upward sides and downward sides of the wings result in lift. In Case 1, the scales of the high-pressure zone and low-pressure zone are much larger than those in Case 2. At the same time, the maximum normalized pressure in Case 1 is approximately 16, while it is only approximately 10 in Case 2. The minimum normalized pressure in Case 1 is approximately −17, while it is only −13 in Case 2. As a result, the lift is greater in Case 1.
For thrust, when the wing is flapping with a greater angle of attack in the down-stroke, the pressure difference contributes to the thrust, which thus largely fluctuates in Case 1, as presented in Fig. 7. A greater α down-stroke leads to a greater pressure difference in the thrust-axis (horizontal direction), so the thrust of Case 1 is much smaller  The vortex fields from τ = 0.38 to τ = 0.48 in Case 1 and Case 2 are shown in Fig. 9. In Case 1, the larger α down-stroke makes it hard for the fluid to flow though the gap between the fore wing and hind wing, so the fluid on the downward side of the wings is compressed. Thus, a high-pressure zone with greater pressure and wider extent is formed on the downward sides of wings, especially near the trailing edge of the fore wing. Then, the fluid is driven forward and flows over the leading edge of the fore wing. Because a higher velocity basically suggests a lower static pressure, a low-pressure zone can be found on the upward sides of wings, especially near the leading edge of the fore wing. In Case 2, the smaller α down-stroke makes it easier for the fluid to flow though the gap between the fore wing and hind wing, so the fluid on the downward sides of the wings is not blocked as much as in Case 1. Thus, the scale of the high-pressure zone on the downward sides of wings is smaller, and the pressure is lower. Meanwhile, the fluid on the downward sides overflows to the upward sides via three routes: 1) through the leading edge of the fore wing, 2) through the gap between the fore wing and hind wing, and 3) through the trailing edge of the hind wing. As a result, the velocity of the upward side is lower in Case 2, and the static pressure is higher in Case 1. The vortex field explains the greater pressure difference across the wings and why the location of the low-pressure zone and high-pressure zone are not the same in the two cases.
In sum, flapping with a greater angle of attack in the down-stroke results in greater lift, while flapping with a lower angle of attack in the down-stroke induces greater thrust but smaller lift. These results are consistent with the experimental results. At the beginning of take-off, the dragonfly requires larger lift, so it will flap the wings with a greater angle of attack in the down-stroke and then switch to a lower angle of attack in the down-stroke mode when it has already left the ground and needs to fly upward and forward.
Effect of phase differences between the forewing and hindwing. Two other cases are also simulated with various phase differences: Case 2 (ϕ = 0 degrees, α down-stroke = 60 degrees) and Case 3 (ϕ = 110 degrees, α down-stroke = 60 degrees). The other conditions are the same.
The lift coefficient (C L ) and drag coefficient (C T ) of wings during one wingbeat, as well as the averaged lift coefficient (C L, avg ) and averaged thrust coefficient (C T, avg ), are shown in Fig. 10. We can speculate that the maximum C L generated by synchronous flapping (Case 3) is larger than that generated by asynchronous flapping (Case 4), and the C L, avg of synchronous flapping is 20.6% larger than that of asynchronous flapping.
In Fig. 10(a), the C L and C T of Case 3 display smaller fluctuations than those of Case 2. The vortex fields in Case 2 and Case 3 from τ = 0.33 to τ = 0.53 are also shown in Fig. 11. In Case 3, both the fore wing and hind wing are flapping in the same direction; as a result, the directions of the aerodynamic forces generated by both the fore wing and hind wing are the same. In Case 4, asynchronous flapping of the fore wing and hind wing leads to opposite directions of the aerodynamic forces generated by the fore wing and hind wing; as a result, the averaged aerodynamic force of both the fore wing and hind wing is more stable.
The thrust generated by synchronous flapping fluctuates with greater amplitude; it is positive when the wings are flapping up and negative when the wings are flapping down. However, in Case 4, the thrust induced by asynchronous flapping almost remains entirely positive. Thus, the averaged thrust is greater in Case 4 than in Case 3.
In sum, the lift generated by synchronous flapping is greater, and when the phase difference between the fore wing and hand wing is larger, a larger thrust is generated. This can also explain why a dragonfly will flap the fore wing and hind wing synchronously at the beginning of take-off, and the phase difference increases when the dragonfly is going to fly upward and forward.
It should be noted that current computational simulations are two-dimensional, whereas three-dimensional aerodynamic effects exist for the flapping of wings. However, this limitation does not affect the results and conclusions drawn above. Comparing the three-dimensional results (Sun and Lan 19 ) with the two-dimensional results (Wang 24 ; Lan and Sun 25 ) of a hovering dragonfly (a situation similar to the current simulations), it can be shown that computed aerodynamic forces show relatively small differences under different dimensionalities. The averaged vertical force coefficient in three-dimensional simulations is approximately 20% less than that in two-dimensional simulations. This discrepancy is consistent with the lift reduction due to the tip vortex, which only exists in three-dimensions. Meanwhile, the time courses of the force coefficients in these simulations are nearly identical. This further indicates that regardless of whether the leading-edge vortex sheds or not, the difference in the vortex structures between two-dimensional simulations and three-dimensional simulations has little influence on the conclusions of this study.
In our simulations, we first focus on the variation in the angle of attack. The alternation of the pitching movement of the wing can be well represented in two dimensions, and the flow phenomenon induced by this alternation is primarily two dimensional. As for the phase difference, Russel and Wang 26 studied its aerodynamic influence on a hovering dragonfly based on two-dimensional flow. The present research employs a similar approach and draws conclusions that are consistent with other two/three-dimensional simulations 19,[25][26][27] . Three-dimensional effects definitely exist in the take-off flight of the dragonfly and will be of interest for future research.

Conclusion
A dragonfly is airborne in approximately one wingbeat and then starts to take-off. The maximum vertical acceleration is close to 20 m/s 2 . The experimental result shows that the acceleration is generated only by aerodynamic force. The phase difference increases from approximately 0 degrees to 110 degrees. The angles of attack in the down-stroke of both wings in the first two wingbeats are significantly greater than those in the later wingbeats, resulting in different levels of aerodynamic lift.
The results of numerical simulation with the 2-D immersed boundary method for various angles of attack in the down-stroke and phase differences between the motion of the fore wing and hind wing show that flapping with a greater angle of attack in the down-stroke generates greater lift but smaller thrust; on the other hand, flapping with a lower angle of attack in the down-stroke induces greater thrust but smaller lift. Synchronous flapping (with zero phase difference) generates greater lift but smaller thrust. With the increase in the phase difference, larger thrust and smaller lift are induced.
Dragonflies are smart and good at using the best method of flapping. At the beginning of take-off, this requires rapidly increasing lift, so the flapping mode is synchronous with a large angle of attack in the down-stroke. After leaving the ground, smaller lift but larger drag is required to move forward, so the mode switches to flapping with increasing phase differences and decreasing angles of attack in the down-stroke.

Methodology
Experimental observation. Preparation of Animals. Pantala flavescens belongs to a widely distributed dragonfly family (Libellulidae) and is considered to be the most widespread dragonfly on the planet. Dragonflies were collected from a pond and lawn in Beihang University between 30 July and 16 August 2015. After this, the dragonflies were stored temporarily in a cool dark storage container and then moved into the indoor laboratory.
In the laboratory, the dragonflies were first put in a refrigerator at a low temperature (~0-2 °C) for 1.5-2 hours. Then, the dragonflies were taken out and kept under anesthesia for approximately 5 minutes at room temperature. In this period, the featured parameters of dragonflies (such as weight and morphological parameters of the wings) were measured. For three-dimensional reconstruction, some black dots were marked on the wings (plotted in Fig. 12). All subjects were tested within 5 hours to ensure their activeness.
Experimental facilities. As shown in Fig. 13, a 78 × 78 × 100 cm 3 observation box made of iron frames and specular glass walls was used to limit the flying space of dragonflies. Two 1000 W photographic lamps were also mounted for illumination compensation in the filming area. Meanwhile, two high-speed cameras (Olympus i-SPEED TR, 1000 frames-1, shutter speed 1 ms, resolution 1280 × 1024 pixels) were used to film the take-offs of dragonflies.
Experimental procedure. Because the objective of this paper is to study the voluntary take-off of dragonflies, the subjects were individually put on the top of a support in the observation box. Then, videos were recorded as soon as each dragonfly moved its wings. For each dragonfly, the process of voluntary take-off was repeatedly captured 2-3 times. The total number of dragonflies tested in this experiment was 30. All the recorded voluntary take-off processes were almost the same, so the results shown in this paper of one dragonfly are representative. Three-dimensional reconstruction. First, the optical axes of the cameras were set to be orthogonal to each other using a circular calibrator. Because the video recorded by each high-speed camera captures a two-dimensional plane, a three-dimensional reconstruction was performed to switch the image in the camera coordinates (u, v) to the world coordinates (X, Y, Z) using the linear coplanar method introduced in refs 28,29 .
The conversion formula is shown as follows: where (u 0 , v 0 ) are the coordinates of the main point of the camera, dx and dy are the size of pixels in the horizontal direction and vertical direction, respectively; f is the effective focal length; k is the distortion coefficient; r is the distance from the point to the optical axis; rotation matrix R and translation matrix T are calculated based on the calibration data; and Z c is the field depth of the cameras.
A validation was carried out using a board with known coordinate information. The distance between two adjacent points on this board was 100 mm, as shown in Fig. 14. For the reconstruction, the maximum error was calculated to be 1.134 mm, corresponding to a maximum relative error of 1.13%. In the real test, the largest scale (the length of wings) was approximately 100 mm, so the maximum error was also estimated to be approximately 1.134 mm.
Kinematic parameter calculation. After the three-dimensional reconstruction, all coordinates of the dragonfly body and the marks on the wings can be obtained in the world coordinate system (X, Y, Z). However, the coordinates in a relative coordinate system (X b , Y b , Z b ), located on the dragonfly, are also needed to obtain the kinematic parameters. Figure 15 shows the method we used to transform the world coordinate system to the relative coordinate system. The body centroid of the dragonfly is defined as the origin point (O b ) of the relative coordinate system located on the dragonfly. Please note that the body centroid of the dragonfly is reasonably defined as the crossing point between the body axial line (from the head of the dragonfly to the tail) and the line connecting the root points of the right forewing and left forewing 18 . After that, the positive direction of the X-axis in the relative coordinate system is defined as the direction pointing toward the tail of the dragonfly and named the X b direction. The Z-axis is set to be perpendicular to the X b axis through the root points of both fore wings. Next, the Y-axis is readily defined because the relative coordinate system is a right-handed system.
To describe the wing kinematics of dragonflies, the method given by Ellington 30 was followed. For the stroke plane of the fore wing, the coordinates (only X b and Y b ) of both fore wing tips in one wingbeat cycle were first marked on the plane of symmetry (plane Z b = 0). Based on these projected points, a line can be determined by the linear regression method. As a result, the stroke plane of the fore wing is defined based on the wing tip line and the root of the wing. For the stroke plane of the hind wing, the same procedure was repeated.
For each stroke plane, the angle between the stroke plane and the body axis is a constant, which is identified as the stroke plane angle β shown in Fig. 16a. In Fig. 16b, a line is drawn between the wing root and wing tip. Let (X, Y, Z) be a reference frame with the origin located at the wing root and the X-Y plane coinciding with the stroke plane. The orientation of the wing with respect to the stroke plane is determined by three Euler angles: the positional angle (ϕ), stroke deviation angle (θ) and pitch angle (ψ), where ϕ is defined as the angle between the Y-axis and the projection onto the stroke plane of the line joining the wing base and wing tip, θ is defined as the angle between the line joining the wing base and the wing tip and its projection onto the stroke plane, and ψ is defined as the angle between the local wing chord and line l (l is perpendicular to the wing span and parallel to the stroke plane). The angle of attack of the wing (α) has the following relationship with ψ: in the down-stroke, α = ψ; in the up-stroke, α = 180 degrees-ψ. were also carried out in this paper. The numerical method is introduced in detail as follows: Immersed Boundary Method. The use of the orthogonal Cartesian coordinates is the most remarkable feature of the immersed boundary method (IBM), which makes the whole grid remain motionless regardless of the movement of the boundary. A source term is introduced into the governing equations to describe the effect of the boundary. The IBM was first applied to simulate cardiac blood flow in 1972 31 . Because of the great simplification of grid generation compared to traditional methods, especially in cases with transforming and moving boundaries, the IBM has been greatly developed and widely applied. Over the past 40 years, researchers have found a better way to address the boundary issues through upgrading and optimization in the IBM [32][33][34] .     where A f and A h are the flapping amplitude of the fore wing and hind wing, respectively; f is the flapping frequency; and ϕ is the phase difference between the fore wing and the hind wing. For the simulations, the coordinates of the center of the model wings can be written as follows: where (x fo , y fo ) and (x ho , y ho ) are the initial positions of the fore wing and hind wing, respectively, and β f and β h are the stroke plane angles of the fore wing and hind wing, respectively. According to the variation in the attack angle observed in real dragonflies in the experiment, the attack angle remains a stable value except at the beginning and the end of both up-strokes and down-strokes.
The rotational movement of the fore wing and the hind wing are described by Equation (6) and Equation (7), respectively, where α is the angle between the long axis of the ellipse and the horizontal plane, α 0 is the initial rotation angle, α m is the stable attack angle in the middle of flapping, and ϕ is the phase difference between the fore and the hind wing.