Kinematic and Aerodynamic Investigation of the Butterﬂy in Forward Free Flight for the Butterﬂy-Inspired Flapping Wing Air Vehicle

: To ensure the stability of ﬂight, the butterﬂy needs to ﬂap its wings and simultaneously move its main body to achieve all kinds of ﬂying motion, such as taking off, hovering, or reverse ﬂight. The high-speed camera is used to record the swing of the abdomen, the movement of the wings, and the pitch angle of the body for butterﬂies during their free ﬂight; the comprehensive biokinetic observations show that the butterﬂy’s wings and body are coupled in various ﬂight states. The swing of the abdomen and the ﬂap of the fore wing affect the pitch motion signiﬁcantly. For theoretical analysis of the butterﬂy ﬂight, a three-dimensional multi-rigid butterﬂy model based on real butterﬂy dimension is established, and the aerodynamic of the butterﬂy ﬂight is simulated and analyzed via computational ﬂuid dynamics methods to obtain an optimal kinematic model of butterﬂy forward ﬂight. Moreover, the formation and development of three-dimensional vortex structures in the forward ﬂight are also presented. The detailed structures of vortices and their dynamic behavior show that the wing’s ﬂap and the abdominal swing play a key role in reorienting and correcting the “clap and peel” mechanism, and the force generation mechanisms are evaluated. The research indicates that longitudinal ﬂight performance is mainly related to the kinematic parameters of the wing and body, and it can lead to the development of butterﬂy-inspired ﬂapping wing air vehicles. of 3D ﬁeld 3D structures of LEV, WTV, and TEV The formation and development of the vortices further to the high-lift mechanism during the free forward ﬂight. A pair of leading-edge vortices is attached to the leading edges of the left and right wings throughout the downstroke. Pressure distribution on butterﬂy wings and body surface reveal unsteady The simulation results of the instantaneous aerodynamic mechanism corresponding


Introduction
The bionic study of insect flight has lasted for a long time, and common characteristics of insect flight have been mentioned in many studies, such as elastic and thin wings and high flapping frequency (10-500 Hz) [1]. Common insects, such as bees, dragonflies, and flies, all conform to these characteristics, while butterflies' wing-flapping frequency is approximately 10 Hz, which is much lower than many other insects [2]. Moreover, the butterfly has a larger wing area and a lower aspect ratio compared with other insects. The flight-related structures of insects are mainly concentrated on the wing-thoracic segment. The muscle groups control the wings to twist and oscillate, thereby generating sufficient lift and thrust. Furthermore, the wings are crisscrossed with wing veins, which support and reinforce the wings similar to a skeleton, as well as can be twisted to change the direction of flight [3]. Tanaka and Shimoyama [4] recognized the effect of the inner veins on the flight performance. Wootton [5] suggested that insects' flying capabilities and various flight skills rely on various wing types and complex wing movement patterns. To adjust the velocity, direction, and attitude of flight, some insects change the movement mode of wings, while others control the flight by changing the relative motion of various parts of themselves.
Take a dragonfly as an example. Its abdominal segment has a certain effect on flight control because its long abdomen can bend or even curl [6]. The butterfly's flight state seems to be unstable, but it can complete the flight target, which is a unique phenomenon in insect flight.
In the 1980s and 1990s, early observations of the flight mechanism of dragonflies, butterflies, and other insects were conducted by Ellington, Dickinson, and Dudley, among others [7][8][9][10][11][12][13][14]. Ellington et al. revealed the effect of leading-edge vortices (LEV) on the high lift of insects [7]. With the development of high-speed photography technology, Wing-motion of hovering small fly or artificial flapping wings have even been measured and flows of the wings calculated numerically [15,16]. It is known that a butterfly flaps downward to generate lift force and upward to generate thrust force. The phenomenon of high lift in insect flight was analyzed qualitatively through experimental fluid methods. "Rotational circulation","wake capture","dynamic stall or the delayed stall", and "Clap and fling" mechanisms were discovered and studied successively [9,[17][18][19]. The kinematics of the wings and body were measured by Sunada et al. [20] via a high-speed video camera during the take-off of a butterfly. They developed a vortex method to compute the aerodynamic force on the flapping wings. The results indicate that the two wings' opening motion generated the first peak of lift-off.
In terms of biomorphological features, the butterfly is mainly composed of the five parts wing, head, thorax, abdomen, and foot, and its flight behavior is relatively complicated. Current research indicates that a butterfly's flight has the following kinematic characteristics: (1) The flapping angle has low frequency and large amplitude during wing flapping. (2) Lead-lag angle: In the downstroke, "lead-lag motion" occurs with the forewings flapping. In other words, the forewings sweep forward along the longitudinal axis of the body to reduce the overlapping area of the fore and hind wings. The results are reversed in upstroke. (3) Feathering angle: The twist angle of the wing is considerable. Furthermore, these characteristics have a certain relationship with the flying ability of the butterfly [21]. At present, the aerodynamic study on butterfly flight is less than that of flies, bees, moths, etc. [2]. Simultaneously, as the wings of the butterfly grow in the thorax's position, the flapping motion of the wings will cause the thorax to produce the pitch oscillation with the same frequency. In most flight states (takeoff, hovering, forward flight, etc.), the abdomen will swing actively with the body's pitch oscillation to stabilize its longitudinal attitude and even influence its maneuvering flight [22,23]. Compared with most birds and insects with a high aspect ratio, the flight behavior of butterflies is more complicated, and the research on its flight mechanism is challenging.
Srygley et al. [24] and Okamoto et al. [25] used alive insects to conduct flight experiments. Bomphrey et al. installed a device that can release smoke on the leading edge of the flying bumblebees' wings and recorded the formation and changes of airflow during the flapping process with the high-speed camera [26]. In these experiments, the constraint of the rope has a great influence on the insect's flight, and the measurement errors of the observation results cannot be ignored. Besides, the measurement has limitations due to living insects' use as experimental subjects in the experiment. For example, it is impossible to measure the force applied to the wings' roots, and it is difficult to confirm the influence of different vibration modes of the wings on the lift. Some scholars have proposed quasi-steady aerodynamic methods to estimate the aerodynamic forces and moments of butterflies flying forward or other insects hovering [20,27]. The lift and thrust mechanism generated by the flapping wing has been analyzed partially, which is influenced by the swinging effect of the butterfly's abdomen during its flight [28][29][30]. Furthermore, it requires more detailed observation and analysis.
Based on the previous studies, our goal is to study the kinematics and aerodynamics of a butterfly in forward flight. The flapping flights of relatively large-size butterflies are chosen to be observed and analyzed. High-speed cameras are arranged to capture the highdefinition forward flight images of butterflies and track the spatial trajectory of the feature points on the butterfly. The kinematics model of the butterfly with wing-body coupling is established to quantify the flight parameters. Based on these results, a three-dimensional multi-rigid butterfly model based on real butterfly dimensions is established. The lift and thrust characteristics of the butterfly are simulated and analyzed by computational fluid dynamics (CFD) methods. Bilateral and symmetrical flapping in forward flight is numerically simulated to clarify the longitudinal flight mechanism of the butterfly. The detailed structure of vortices and their dynamic behavior are elucidated, and their relations to force generation mechanisms are evaluated.

Materials and Equipment
To study the special flying skills of the butterfly, it is necessary to conduct detailed observation and quantitative analysis of the free flying state of the live butterfly. Chilasa clytia were collected from a breeding base in southern China. All adult butterflies are fed with fresh, rotten apples in the laboratory and kept active in a container with the temperature ranging from 25 to 32 • C and the humidity ranging from 40% to 70%.
Morphological parameters of the butterfly (Chilasa clytia) are summarized in Table 1, including the total mass (m), the thoracic mass (m t ), the abdominal mass (m a ), the thoracic length (L t ), the thoracic width (w t ), the abdominal length (L a ), the abdominal width (w a ), the wings mass (m w ), the wing-tip length (L tip ), the wing area (2S w ), the wing loading (mg/(2S w )), and the chord length of the wing (c). The mass of the butterfly is approximately 0.35 g. The wingspan (S) is about 8.29 cm, and its chord length is about 3.55 cm. The aspect ratio is expressed as AR = S/c. AR of the butterfly is about 2.3. The biomorphological characteristics of the butterfly are obviously different from other insects such as cockroaches, hawk moths, and bee flies.
Because the total mass of the live butterfly can be measured directly, the thorax masses and the wings subtracted from the total mass give the abdominal mass. Note that the actual value of the abdominal mass largely varies during the observation. The mass distribution of various body parts of the experimental butterfly is calculated. m t is similar to m a , and they account for about 30% of the total mass, respectively. m w accounts for about 15% of the total mass. Thus, it can be concluded that butterfly wings have very small surface mass density.

Parameters
Value Compared with other flying insects, the butterfly has a larger physical size (the maximum span is more than 10 cm). Large flying space is required to ensure sufficient time to record the flight, so the transparent cube container made of grass with dimensions of 0.4 m × 0.4 m × 0.4 m is adopted. The flight velocity of the butterfly is high, and the frequency range of flapping motion is generally about 9-12 Hz. Therefore, enough light sources are arranged to increase the quantity of light entering the image sensor to ensure the sharpness of the image. The observation experiment system is shown in Figure 1d. The red dot is the initial release position of the butterfly, which is 30 cm high from the ground. Horizontal body plane  (d) experiment system; and (e) three coordinate systems to describe the motion of a butterfly, namely the inertial coordinate system (X, Y, Z), the body-fixed coordinate system (X b , Y b , Z b ) and wing-fixed coordinate system (X w , Y w , Z w ).
Two kinds of light sources are used in biokinetic observation. Light A is a QL-1000 LED continuum light source, and the output power is 1000 W, which is arranged in front of the cube container. Light B is OPT-FLP 125 mm × 103 mm surface light source, and the output power is 24 V/6 W, which is arranged directly above the cube container, to make sure the subject will not be shaded as much as possible. To observe the kinematic characteristics of various parts (abdomen, thorax, wing, etc.) in various flight states more clearly, the high-speed camera Motion BLITZ Cube4, which is made in Germany, is used for the biokinetic observation of butterfly. The frame rate is up to 93,000 fps (frames per second), and its output is black and white video. The resolution range of the captured image is 1280 (H) × 10 (V)-1280 (H) × 1024 (V) pixels. The shutter speed can be set to 2-1000 µs.

Coordinate System of the Flight Kinematics
In this study, the Cartesian coordinate system is introduced to describe the wing and body kinematics of a butterfly. Three coordinate systems are applied to define the axis, the plane, and the characteristic angles related to the butterfly motion in Figure 1e: The first is the inertial coordinate system (X, Y, Z); its origin (O) is set at the initial center of gravity (COG), which coincides with initial release position. The X-axis is defined therein as the direction in which a butterfly moves forward, while the Y-axis is the direction opposite to that of gravity. The thorax rotation around the Z-axis is described with pitch motion, which is defined as the body pitch angle between the centerline of the thorax and horizontal plane (X-Z plane). For instance, if the pitch angle is zero, the body will lie in horizontal plane; when it is 90 • , the body will be vertical.
The body-fixed coordinate system (X b , Y b , Z b ) is adopted to describe the body motion. The origin of the body-fixed coordinate system (O b ) is set at the COG of the butterfly, which translates and rotates with the body. During the forward flight, the flight velocity of the butterfly is expressed as V, and the butterfly's displacement is decomposable into two directions in the inertial coordinate system, the Xand Y-axis. The velocities in the X and Y directions are horizontal velocity v x and vertical velocity v y , respectively. The symmetry axis of the butterfly thorax is defined as the thorax axis, which coincides with X b -axis. The symmetry axis of the butterfly abdomen is defined as the abdomen axis. X b -Z b and X b -Y b denote the horizontal and vertical body planes, respectively. The abdomen rotates around the Z b -axis, and the angle is defined as abdominal swing angle.
The wing-fixed coordinate system (X w , Y w , Z w ) is adopted to describe the wingflapping motion. The motion of the wings is connected and simultaneous with the butterfly body in real flight. The origin of the wing-fixed coordinate system (O w ) is set at the wing root and translates and rotates with the body. The butterfly flaps its wings around X w -axis, which is simplified to be parallel to the thoracic axis in this paper; the flapping angle is 90 • when the two wings clap together.
The plane swept by the connection between the wing tip and the wing root during the fight is defined as the stroke plane of the butterfly wing. The angle between the stroke plane and the horizontal X-Z plane is defined as β. In this study, the thoracic pitch angle, the abdominal swing angle, and the approximate flapping angle are denoted by θ p , θ a , and θ f , respectively. This paper mainly considers the longitudinal flight of butterflies, so the roll and yaw motions are ignored.

Movement Tracking of Butterfly Flight
We characterized and analyzed the motion and postures of the experimental butterfly flying horizontally in the cube container. The forward flight might be defined as translation horizontally in space. Because of the difficulty of recording the specific forward flight of a butterfly in the laboratory experiment, we redefined the forward flight such that the ratio between the vertical and horizontal translation of the COG is smaller than 0.15, |dy/dx| < 0.15, within a cycle. The wings and body movements were recorded with a high-speed camera and operated at the frame rate of 1000 fps with pixel size 1280 × 1024. To reduce the calculation amount of image processing for flight trajectory, the collected video is framed with pixel size 500 × 900 so that the calculation amount is reduced by about 50% and the processing time is reduced by more than half. The image-based feature point recognition method is utilized to track and analyze the motion.
The four feature points for the observation records of the spatial displacement are shown as purple thoracic COM (x t , y t ), white joint point (x j , y j ), blue abdominal end (x a , y a ), and red wingtip (x wt , y wt ) in Figure 1c. The thorax axis and the abdomen axis are shown as solid green and light blue lines, respectively. Since the wing root is very close to the thoracic COM, the wing root and the thoracic COM are considered to coincide during the measurement. The joint angles of the thorax and abdomen during the free flight are observed. θ p and θ a are obtained from the geometric relationship, and numerical values of the feature points coordinates are derived as in Equations (1) and (2) The finite difference algorithm adapts to calculate velocity and acceleration in the X-axis directions (the same with the Y-axis directions) is given by where v x i and a x i are the velocity and acceleration of the feature points in the X direction at frame i, respectively. x i is the position coordinate value of the feature point on the ith frame image. dt is the time interval between adjacent frames of the video. Then, the derived velocity and acceleration are low-pass filtered to remove noise. In the experiment, the instantaneous wingtip velocity was calculated as follows: and v wt y represent the velocity of the wing tip in the X and Y direction, respectively.

3D Multi-Body Model of a Butterfly
There are a fore wing and a hind wing on each side of the butterfly body. To further simplify the model, the fore wing and hind wing are modeled as a unified wing. Figure 1a shows that the experimental butterflies (Chilasa clytia) consist of four parts, namely a spheroidal thorax, a spheroidal abdomen, and one pair of thin flat wings, and they are set as rigid bodies in the simulation. The thorax and abdomen are connected by joints with one degree of freedom (DOF). The thorax and left wing are connected by a joint with three DOFs, as are the thorax and right wing. The parameters of the real butterfly are used as the dimensions of the simulation model for CFD.
The wingspan and chord length of the simulation model is set to 8.5 and 3.6 cm. The positions of the leading edge and the trailing edge on the wings are shown in Figure 1b. Based on the two-dimensional images of a real butterfly, a simplified butterfly model is created. All parameters defining butterfly simulation model geometry are provided in Table 2. Table 2. Parameters defining the simulation model of the butterfly.

Computational Setup for CFD Simulation
To clarify the influence of butterfly kinematics parameters on free flight, the CFD simulations are conducted in this paper. The direction of the mainstream is taken as the positive X direction, the vertically upward direction is taken as the positive Y direction, and the horizontal direction is taken as the Z direction. The directions are expressed by the Cartesian unit vectors e e e x , e e e y , and e e e z . The equations governing the flow field include an equation for the three-dimensional, transient, and incompressible continuity and the Navier-Stokes equation, ∇ · u u u = 0 (5) where u u u is the velocity vector of the flow in the Cartesian coordinate system, defined as u u u = u 0 e e e x . p is the pressure. The fluid is considered to be 20 • C. The values of the air properties are adopted from those at 20 • C: the density of air ρ air = 1.205 kg/m 3 and the kinematic viscosity υ air = 1.5 × 10 −5 m 2 /s. It should be noted that the gravitational term does not appear in Equation (6). The velocity is constant on the inflow boundary plane: V ∞ , incoming flow velocity. The streamwise length of the computational domain is L x = 0.8m, which is approximately 20 times as long as the aerodynamic chord length. The boundary conditions of the flows are shown in Figure 2a. The periodic boundary condition is employed for the transverse boundaries, and, on the other sides, parallel to the longitudinal X-axis, the no-slip condition is used. The transverse length of the computational domain is L y = 0.4 m, L z = 0.5 m, which is approximately five times as long as the wing-tip length. The results are not affected by the values of L x , L y , and L z . The average wing-tip velocityV tip equals the reference velocity defined as V re f = 2A f πr f /180 • . A f is the stroke amplitude, with the value about 120 ± 10 • , f is the flapping frequency, r is the single wing length which equals L tip /2. The governing parameter of the above equations is the Reynolds number (Re) given by where L re f is the reference length which is equal toc.c is the mid-chord length (2.5 × 10 2 m).
It should be noted that the definition of the Re is the same as that used in [2,31].
In this simulation case, the Re during butterfly forward flight is about 3000. Two key parameters of the flapping wing propulsion are the Strouhal number [32] and the advanced ratio, which are defined as  The aerodynamics model of the butterfly is simulated numerically with the Lattice Boltzmann method (LBM) [33], which is applied with a commercial LBM code: Next Limit Technologies' XFlow. Unlike traditional CFD software, XFlow provides a proprietary fully Lagrangian particle-based kinetic solver, which dramatically reduces the time to calculate for a given level of solution accuracy. The LBM proves to be efficient and reliable for aerodynamic force, torque, and vortex wake structure simulation under various conditions [34][35][36]. The kinematic law of the butterfly is adapted to the multi-body model by the kinematic simulation software (ADAMS). Firstly, the three-dimensional model of the butterfly is established in the modeling software Solidworks. Then, the model is imported into ADAMS. In ADAMS, kinematics rules for the rigid body parts are supplemented and then the ".res" file is generated. Finally, the three-dimensional model and motion rule file of the butterfly are imported into XFlow. The fluid conditions can be set and then carry out simulation data analysis.
The spatial resolution of the far field (Resolved scale in XFlow) is set to be the largest level of the lattice in adaptive refinement mode. The smaller the resolved scale is, the more accurate the numerical calculation is. The simulation time also increases with the number of elements being calculated. The resolution of the near wall is set to smallest to refine dynamically moving geometries. The dynamic domain topology will perform octree stratification according to the several resolution settings to get the ideal computational lattice, as shown in Figure 2b. The resolution of the wake is set to the same value as the resolution of the near wall for getting the clear wake flow.
In this study, the time step of each flight case is set to 0.001 s, and the resolution scheme is carefully selected, as shown in Table 3. Note that "s" represents the initially tested resolution. Figure 2c,d shows the comparison of aerodynamic forces during the four flapping strokes from three resolution set-ups (coarse, medium, and fine). The differences of the average values of the horizontal and vertical aerodynamic forces between the medium-resolution (adopted in this paper) and the fine-resolution are about 2% and 10%, respectively. In XFlow, it is necessary to ensure that the stability parameter of the simulation is low, and the relatively fine domain structure can be obtained. Considering the calculation time and accuracy, the medium-resolution scheme of s/2 is chosen as the simulation condition. The flight motion of the butterfly mainly couples the motions of the wing flapping and the thorax pitching. The abdominal motion, which actively changes the pitch Angle of the thorax, is important for the stability of the flight. By pitching the body, the butterfly manipulates the stroke plane angle to adjust the direction of the flapping motion. The movement of wings is accompanied by clapping and peeling during the downstroke (DS) and upstroke (US), respectively.
The flapping motion and the pitch motion of the real butterfly and the simulation model change periodically, as shown in Figure 3. The swing of the abdomen directly affects the pitch angle of the butterfly thorax. Thus, the wing kinematics can be expressed as the flapping motion θ f (t), the pitch motion θ p (t), and the phase difference γ between θ f and θ p . The wing and body kinematics of four consecutive flapping cycles were recreated with our programs. To input the flight kinematics model of butterflies recorded from observation into a simulation, three cosine wave functions served to approach the flapping motion and pitch motion the same as in the experiment. The functions are expressed as in which A f , A p , and A a are the amplitudes of θ f (t), θ p (t), and θ a (t) from the experiment. The aerodynamic forces are calculated with mathematical sets of the kinematic parameter

Wing and Body Kinematics of Butterfly Forward Flight
Our aim in this paper is to present the clearest forward flight and reveal the special flying skills of the butterfly (Chilasa clytia). Ten videos of the butterfly free flying about 300 ms long were filmed. Two of them were selected as the two cases (denoted as VC1 and VC2) of forward flight because the displacement of the butterfly centroid in the X-axis and Y-axis directions satisfied |dy/dx| < 0.15.
For the two cases, the analyzed kinematics results of the butterfly are shown in Table 4, including the total time of forward flight (t total ), the average stroke amplitude (Ā f ), the average flapping frequency (f ), the average pitch angle (θ p ), the maximum pitch angle (θ p max ), the minimum pitch angle (θ p min ), the average abdominal swing angle (θ a ), the maximum abdominal swing angle (θ a max ), and the minimum abdominal swing angle (θ a min ). Two complete flapping cycles were intercepted from the two videos, named Stroke #1 and Stroke #2, which correspond to the subscript numbers, respectively.  ID θ 1p max θ 1p minθ1 p θ 1a max θ 1a minθ1 a θ 2p max θ 2p minθ2 p θ 2a max θ 2a minθ2 a t totalĀ  For each of the flights, the recordings of 2-3 wing flapping cycles were digitized. Snapshots of a butterfly in forward flight VC1 and VC2 from the front-view camera are shown in Figure 4a,b, respectively. Because the displacement along the Z-axis is very close to zero throughout the entire process of forward flight, only the displacements along the X-axis and Y-axis are plotted. The coordinates at 0 ms are set to (0, 0, 0).  After tracking all the feature points based on all recorded frames, the displacements and kinematic curves of the butterfly in the relative coordinate system were obtained from VC1 and VC2 and plotted in Figures 5 and 6, respectively. All the stroke moments were also labeled. Grey shading indicates the DS, while white shading denotes the US. The purple shading indicates the gliding flight from the end of the US to the beginning of the next flapping cycle. After being released from the initial release position, the butterfly first flaps its wings slightly for one or two strokes with a small reduction in flight altitude, and its average flight speed increases accordingly; then, the average pitch angle increases, which is equivalent to an increase in the angle of attack of the body. The butterfly starts to climb, and its average flight speed decreases accordingly.
The vertical acceleration (a y ) reaches the maximum value at the mid-downstroke of the DS, at which time a lift peak appears, and almost no lift is generated during the US. For almost the entire downstroke, a y is positive because the instantaneous vertical aerodynamic force of the butterfly is lift, while the horizontal acceleration (a x ) is positive because the instantaneous horizontal aerodynamic force is drag. a x reaches the maximum value at near the mid-upstroke of the US, at which time a peak of the thrust appears. The kinematics analysis result shows that the phase difference between butterfly flapping motion and pitching motion is close to π/2. Under different flight conditions, the phase difference between the pitch angle and the swing angle varies greatly, which is about π in our observation.  During the flight, the glide flight mode is used to reduce flight energy consumption. This mode often appears between the end of the upper flap and the beginning of the next flapping cycle. At this time, the wings are almost closed on the back, the projection area in the flight direction is the smallest, and the air resistance is also small. Therefore, the main body has a certain height drop during the gliding process due to its own inertia, and the flight velocity will increase accordingly. Due to the unstable flight of the butterfly, it is difficult to capture the state of the butterfly's continuous forward flight under free flight conditions. In this experiment, the butterfly is released from a fixed position to guarantee that the butterfly's head faces the −X direction of the space.
For VC1 and VC2, the flying behaviors of the experimental butterfly are similar. After the butterfly is released, it accelerates to fly forward in the initial stage, and the velocity and flight trajectory have fluctuation characteristics with the same frequency as the flapping motion. When the flight speed increases, the amplitude of the abdominal swing and the average swing angle increase, with the increase of the average pitch angle of the body. Obviously, the upper and lower limit angles of the pitch movement both increase, as shown in Figures 5d and 6d. At the end of the videos, the butterfly has completed the efficient climb flight. In summary, the flight of a real butterfly is complex and unstable, and detailed measurement and analysis of its kinematic parameters under different flight conditions are very meaningful. In the following, the CFD method is used to quantitatively analyze several key kinematic parameters that affect butterfly flight performance.

Effect of Kinematic Parameters and Computational Fluid Dynamics Simulations
In the simulation, the flapping frequency of the butterfly is set to 10 Hz. The aerodynamic force generated by flapping wing at a certain inflow velocity can be obtained through simulation analysis. It displays the overall forces (F x , F y , F z ), in global axes, exerted by the fluid on the entire geometry model and the corresponding force coefficients (C x , C y , C z ): where F x and F y are the instantaneous horizontal and vertical aerodynamic force in downstroke and upstroke, respectively. When F x is negative, it indicates thrust; when F x is positive, it indicates resistance.F x andF y are defined as the average horizontal and vertical aerodynamic force (F aero = 1 T t+T t F aero (t)dt), respectively. The vertical force is the lift force, while the horizontal aerodynamic force is the the resultant force of thrust and resistance. A re f is the reference area which equals the projected area of the wing. S w is the area of one wing.
The kinematic parameters in Equation (8) would change regularly to reveal the effect of the kinematics model on butterfly flight. Different pitch motions may affect the flying velocity of the butterfly, while the fixed flying velocity is considered here. Based on previous observations and combined with empirical analysis, the initial kinematic simulation model can be obtained as Equation (9). Four kinematic parameters including B 0 , V ∞ , γ, and A f with greater influence on butterfly flight were selected as research content. To verify the ideas proposed in the experiment, the following four cases are simulated and analyzed. All simulations start flapping downward from the upper limit position. The four cases of set-ups are shown in Table 5. First, Case 1 is simulated as the free flight without rotation. Otherwise, attitude control is assumed to be perfectly achieved and does not affect the translational motion of the body. The incoming flow velocity is set to 1 m/s according to the actual flying velocity of the butterfly in the horizontal direction.
In this case, the thoracic pitch angle is a fixed value (θ p (t) = B 0 ) to verify that body pitch is necessary for the flight. Moreover, this means that the body of the model does not rotate in the pitching direction and the yawing and rolling directions. B 0 is taken as 4 • , 14 • , 24 • , 34 • , 44 • . The simulation results are as follows (Figure 8). Under the condition that the pitch angle value is small (θ p = 4 • , 14 • ), the positive and negative peak values of the lift generated by flapping are basically equal in magnitude, and the drag is relatively small. As the pitch angle increases within a certain range, both the positive peak value and the mean value of the lift increase. In comparison, the resistance increases during this period because of the increase in the average windward projection area of the butterfly. In summary, the case with a fixed pitch angle (without wings rotation) shows that the maximum average lift of the butterfly is less than 60% of the butterfly's gravity. Thus, the pitching motion of the body is extremely necessary for butterfly flight.

Effect of the Flight Velocity in Forward Free Flight
This simulation aims to investigate the effect of kinematics parameters of the model with rotation in addition to the translational motion. In Case 2, V ∞ is taken as 0.6, 0.8, 1.0, and 1.2 m/s. A kinematic parametric study on flight velocity is presented to explore its influence on aerodynamic performance. We simulate free flights of the butterfly model with rotation (body pitch motion). The kinematics of thoracic motion is simplified to a time-varying cosine function. The initial position and amplitude of the thoracic pitch angle are certain values. The simulation model of the body can rotate in the longitudinal direction, while it does not rotate in the yawing and rolling directions. The simulation results are shown in Figure 9.  As the velocity increases, the lift and drag generated by flapping wing increase synchronously, and their instantaneous peaks increase accordingly.F x andF y of the butterfly at the velocity of 1.0 m/s are 4.7 × 10 −4 and 2.7 × 10 −3 N (approximately 80% of butterfly gravity), respectively. The lift reaches the instantaneous maximum 9.5 × 10 −3 N at the mid-downstroke of the DS (approximately 0.025-0.03 s). The thrust reaches the instantaneous maximum 1.7 × 10 −3 N at the mid-downstroke of the DS (approximately 0.075-0.08 s). The effect of flight velocity on lift and drag is relatively linear within the test velocity range. According to previous kinematic observations, the butterfly has fluctuations in the Y-axis direction, and the velocity in the X-axis direction is also time-varying. v x and v y will be set as a periodic change function in further research.

Effect of the Phase Difference between Flapping and Thoracic Pitch Motion
The coupling of the butterfly body and wings is considered in Case 3. γ is taken as π/6, π/4, π/3, π/2, and 2π/3. The simulation results in Figure 10 show that the lift generated by flapping wing always increases, but the drag decreases first and then increases as γ increases. when γ ≥ π/2, the generation of the drag peak is postponed backward, and the generation time is shifted by about 0.02 s compared to the smaller γ. To achieve forward flight, the butterfly needs to have a certain thrust (low resistance) to get sufficient lift. That is, it is necessary to choose a suitable kinematic parameter γ.

Effect of the Flapping Angle in Forward Free Flight
In Case 4, different wing-flapping motion amplitudes are shown: the flapping balance position is determined, while the upper limit and the lower limit positions during flapping are different. The incoming flow velocity, in this case, is set to 1 m/s. The kinematics rules are set by ADAMS.   Figure 11 shows that the appearance of the peak of the lift and thrust is advanced in time, andF x of the butterfly gradually decreases as the flapping angle increases. That is, the resistance decreases, and thrust increases. It is worth mentioning that the lift force increases with the increase of the A f within a small amplitude range, but, when it is increased to larger than a certain angle (A f = 60 • , in this study), the lift force suddenly decreases to less than 50% of gravity. Thus, it is not that the larger is the A f , the greater is the average lift. The result reveals that the butterfly has a relatively large average lift (close to the butterfly's gravity) when A f is 60 • , consistent with the observation results. The butterfly needs to determine the amplitude of each flapping according to the different lift requirements.
A v e r a g e A e r o d y n a m i c F o r c e ( N )

Vortex Formation and Development
The butterfly uses the chest muscles to drive the wings to stir in the viscous air. The air around the wings generates vortices, which in turn generate enough lift to overcome its own gravity. When the wings flap downward, the body floats upward. When the wings flap upward, it mainly generates thrust to fly forward. The formation and development of 3D vortex structures in one flapping cycle of forward flight are presented in Figure 12. Leading-edge vortex (LEV) and wing-tip vortex (WTV) appear on the dorsal side of the wings, and they are connected at the corner of the wings. The trailing edge vortex (TEV) is also clear during the downstroke. The two kinds of LEVs are reported in hawk moth's hovering flight [37]. As a real butterfly, WTV and LEV are a major cause of the lift and thrust generation of the model. During the DS, the airflow flows upward from the leading edge to the dorsal side of the wings, forming a tube-type vortex over the entire leading edge of the wings (from the wing root to the wing tip) [38]. This vortex is a conical spiral along the direction of the leading edge. This spiral vortex is small at the wing root but larger at the wingtip.   [39], as shown in Figure 13b. During the DS, the wings of the butterfly will produce a slanting downward jet. The reactive force of the jet against the butterfly can overcome the butterfly's own gravity. The wake formed by every flapping stroke is mushroom-shaped from the top view (Figure 13a).
At the same time, during the DS, the entire vortex remains parallel to the leading edge of the wings, and is close to the leading edge without detaching (Vortex shedding). From the wing root to the wing tip, the Leading-edge vortex during the DS (LEVd) gradually weakens, disappears, and then is replaced by another newly generated vortex named wing-tip vortex during the DS (WTVd). After the vortex moves to the tip of the wing, it is twisted with the airflow at the wing tip from the ventral surface of the wing and then extends to the tail of the butterfly, so that the cyclones generated by the flapping motion of the two wings form a vortex ring.
Every time the butterfly flaps its wings downward, a vortex ring forms at the leading edge of the wings, and then it rushes out of the cyclone during the US. The separation on the leading edges makes complex secondary flow vortices named secondary flow vortex (SFV) below the LEVd [39]. At the start of the US, the small vortices form and attach near the thorax and wing root. At the end of the upstroke, a relatively small vortex ring forms on the ventral side of the butterfly wings near the leading edge and the wing tip, and the thrust peak appears simultaneously. Figure 14 shows the exquisite 3D vortex structures in the butterfly forward flight. Each DS, a strong oblique downstroke jet is generated, indicating the generation of the lift. A relatively weak oblique upstroke jet will be generated when it flaps the wings upward, and this indicates the generation of the thrust, as shown in Flight experiments have proven that this vortex created a low-pressure area over the wings of the insect, such as cicada in [41]. Under the impetus of the surrounding high pressure, the wings in the center of the low-pressure area will be attracted to rise. Each low-pressure area will be held up by the surrounding high-pressure area so that insects in the low-pressure area will get a lift.
Previously, aerodynamicists suspected that leading-edge vortices generated by insects' vibrating wings could get enough rising power because they did not see the details of the vortices and were conservative. The researchers believe that the significance of this discovery is that, for the first time, the vortices at the leading edge of the insect's wings have been spiraled, conically advancing from the wing root to the wing tip. This movement keeps the entire cyclone stable, delays the separation of the vortices from the wings, and enables the vortices to provide insects with ascending thrust for as long as possible.  Figure 15 shows that the low-pressure regions appear on the dorsal side of the wing, whereas the high-pressure regions appear on the ventral side in the DS. Due to this pressure distribution, the lift force has a large positive value currently. This pressure distribution can be caused by the vortices (WTV and LEV) on the dorsal surface of the wing and the increase in the pressure of the air on the ventral side. The low-pressure region is attached near the leading edge on the dorsal surface of the wings and gradually comes off as it gets closer to the wing tip. It is likely to be created by the LEV. The low-pressure region appears near the wing tip on the dorsal surface of the wing. It seems to be induced by the WTV. On the other hand, the high-pressure region appears on the ventral side of the wing for all cross sections because the downward velocity of the wings is the maximum around this time, and it is likely to be larger than the downward velocity of the ambient air. Therefore, the air in the ventral side is pressed by the wings, and, consequently, the pressure increases in this region. The pressure difference between the upper and lower surfaces of the hind wing relative to the front wing during the entire flight is small, so the contribution of the hind wing to lift and thrust is limited. During the entire flapping process, the pressure difference of the hind wings is always small due to the lower motion velocity of the hind wings relative to the air and the less influence of the leading-edge vortex. Therefore, the fore wings contribute mainly to the generation of lift and thrust, and the hind wings may play a positive role in flight stability.

Discussion
To reveal the unsteady aerodynamics of the butterfly in free flight, kinematic parameters under different flight conditions were obtained through biological observations, and corresponding kinematic curves were fitted. Based on these, the influence of butterfly kinematics parameters on flight performance was studied. The simulation analysis clarifies that the pitching motion of the body is extremely necessary for butterfly flight, and the balanced position of the pitching motion needs to be an appropriate value, which is about 34 • . The effect of flight velocity on aerodynamic forces is relatively linear within the test velocity range. The average lift and drag generated by flapping wing increase synchronously with flight velocity, and the peak values of the instantaneous aerodynamic forces increase accordingly. Besides, as it increases, the lift always increases, but the drag decreases first and then increases. It is necessary to choose an appropriate kinematic parameter to ensure that the lift and thrust are sufficient for flight.
The amplitude of the flapping angle is not as large as possible. The simulation of this research indicates that the butterfly has a large average lift (close to the butterfly's own gravity) and produces a certain thrust when it is 60 • , which is consistent with the observation results. The butterfly needs to adjust the related parameters in real-time in cooperation with the body pitch movement according to different lift requirements for achieving the best kinematic performance.
A detailed computational fluid dynamics simulation analysis was performed by the Lattice Boltzmann method to obtain the kinematic parameters which influence the generation and change of lift and thrust. The optimal kinematic parameters of the butterfly under the determined flight velocity were found. Visualization of 3D vortex structures of the flow field around the model was applied, and 3D structures of LEV, WTV, and TEV were analyzed in this study. The formation and development of the vortices were observed further to explain the high-lift mechanism during the free forward flight. A pair of conical leading-edge vortices is attached to the leading edges of the left and right wings throughout the downstroke. Pressure distribution on butterfly wings and body surface reveal unsteady aerodynamics. The simulation results of the instantaneous aerodynamic curve can confirm the generation mechanism of the corresponding aerodynamic force.
Although the experimental cases are not rich enough, the qualitative analysis of the kinematic parameters of the butterfly forward flight is obtained, and the variation trend of the kinematic parameters affecting the flight mechanism is clarified. More experiments and simulation work of butterfly free flight will be carried out in the future.

Conclusions
Observations by the high-speed camera system show that the kinematics of the butterfly is particularly different from other small insects or very small insects: the AR of the butterflies is lower (even less than 2), the flapping frequency is lower (about 10 Hz), and the wings are not flipped during flapping. The butterfly's body moves up and down with the obvious pitch motion. The high coupling of the body and the wing motion is the key to the lift generation and attitude adjustment during forward flight. Although the body pitch caused by unpredictable disturbances in the air can be compensated by the complex motions of wings such as forward sweep, backward sweep, and torsion, the stability of the longitudinal flight of the butterfly is still reduced when the abdomen is fixed relative to the thorax. Besides, the butterfly glides forward at a smaller open angle with wings close on its back before the DS. Although the glide time is short, it is considered a clever flight strategy to reduce energy consumption.
Through the co-simulation of ADAMS and XFlow, the numerical simulation calculation of the aerodynamic force of the butterfly in the vertical and horizontal directions is realized. Based on the simulation analysis of four key kinematic parameters, the influence of each kinematic parameter on the flying performance of butterflies was obtained. Furthermore, an optimized multi-rigid body flight kinematics model of the butterfly was established. Simultaneously, Lattice Boltzmann method was used to show the vortices structure generated by the butterfly during forward flight. By analyzing the formation and development of vortices, the characteristics of butterfly flight were revealed.
In the future, a flapping wing air vehicle inspired by the flight mechanism of the butterfly will be designed and developed to refer to the similar kinematics equation in this study so that the controllable longitudinal stability of the robot will be achieved.

Informed Consent Statement: Not applicable.
Acknowledgments: The authors thank Longfei Zhao for his help with the numerical method proposed and simulations conducted.

Conflicts of Interest:
The authors declare no conflict of interest.