The Flight Mechanism of a Bird-like Flapping Wing Robot at a Low Reynolds Number

The ﬂight mechanism of a bird-like ﬂapping wing robot at a low Reynolds number was studied in this study for improving the robot performances. Both the physical model and the kinematic model were ﬁrst established. The dynamic model of the robot at a low Reynolds number was built with the RANS (Reynolds-averaged Navier-Stokes) equations and the Spalart-Allmaras turbulence model. The ﬂight experiments were carried out and the results were discussed. Lift and drag coeﬃcient curves show that it generates upward lift and forward thrust in the phase that the wing ﬂaps downwards, the rate of the coeﬃcient curves is the biggest when the ﬂapping direction changes. Pressure contours indicate that small vortexes with high pressure values appear at the wing edges. There are four velocity vortex groups in total at the front and back of the wing in the velocity contours. Some methods for improving the robot ﬂight eﬃciency and the robot strength as well as the stitching position of the robot skin have been obtained from the above results. The methods provide the important guidance for the stable ﬂights of the ﬂapping wing robot with the high eﬃciency.


Introduction
A bird-like flapping wing robot has many advantages compared with a fixed wing robot and a rotary wing robot. Due to "V" formation and some versatile flight modes, it is possible for bird-like robots to fly in an energy-saving way [1][2][3]. A bird-like flapping wing robot has been widely applied in aerial photography, disaster rescue, military operation and the other fields. For example, due to similarities between a bird-like flapping wing robot and a real bird, it becomes possible for the robot Robird to chase away birds from runways by mimicking a predator [4]. Due to its extensive applications, both models and experiments of the robot have been widely studied [5][6][7][8][9][10]. In spite of the fact that some studies have been conducted in recent decades, there are still some challenges from the air flow fields around the robot, which are vital to reveal the flight mechanism of the robot [6,10,11].
Researchers have extensively studied the mechanical models of bird-like flapping wing robots. Bunget et al., designed a bat-like flapping wing robot which was composed of platform body and bat-inspired wings as well as SMA (Shape Memory Alloy) wires, the wings were connected by flexible joints [12]. As shown in Figure 1(a), Chirarattananon et al., designed a flapping wing robot that was equipped with two piezoelectric actuators, and the rotational motion of the wing was realized by a linear actuator and a flexure-based four bar transmission [11]. As shown in Figure 1(b), Leys et al., designed a flapping wing robot, the wing motion was imitated by an artificial wing, which can generate harmonic strokes with a large amplitude and a high frequency [13]. As shown in Figure 1(c), Grand et al., designed a flapping wing robot whose wing can realize an arbitrary trajectory with the high efficiency [14]. In spite of the fact that some flapping wing robots have been studied, there are less attention on the high flight efficiency and the robot strength as well as the stitching position of the robot skin, and it is critical to design the mechanical model of the robot based on the above aspects [15].
Researchers have established the kinematic model of bird-like flapping wing robots in recent decades. Bunget et al., analyzed the joint data which were mainly from image processing and inverse kinematics, and established the kinematic model of a bat flapping wing [12]. In [14], Grand et al., obtained the kinematic data from the wing movements to characterize the kinematic model of a robot. In [16], three set kinematic equations were used to describe the motion of a flapping wing in the primary and secondary as well as third feathers. Dadashi et al., established the kinematic model of a bat-like robot in flight, and represented the kinematics of an articulated flapping wing with the Denavit-Hartenberg convention method [17].
For decades, many researchers have established the dynamic model of bird-like flapping wing robots. In [18], the averaged dynamic model of a flapping wing robot was established with the force-moment method for describing the robot motion. In [19], the dynamic model of a bird-like robot was established with the strip theory to calculate the aerodynamic force in the flight. Based on the Euler method, Chirarattananon et al., established the dynamic model of a robot to obtain the attitude dynamics, and calculated the total torques imposed on the robot joints [20]. Teoh et al., considered aerodynamic dampers to establish the improved dynamic model of a RoboBee robot [21]. In [22,23], some other methods were also used to establish the dynamic model of a flapping wing robot. e above dynamic models are mainly applied to the rigid bodies. It is difficult to analyze the complex air flow fields during the robot flight with the above models. erefore, it is necessary to establish a dynamic model that is closer to the real flight for obtaining the robot flight mechanism.
e simulation experiments of flapping wing robots have been performed with the CFD (computational fluid dynamics) method. In [24], air flow fields of a flapping wing ornithopter were simulated with the CFD method, and the simulation results were consistent with experiments. Su et al., divided the meshes of a bird-like robot and its air flow fields, and used the CFD method to analyze the ground effect on the robot [16]. Moelyadi et al., obtained the velocity contours of a bird-like robot with the CFD method for analyzing the robot motion characteristics [25]. In [26], air flow fields of a robot were simulated with the CFD method, and the aerodynamic performances at a low Reynolds number were obtained. In [27,28], the transient flow fields of flapping wings were analyzed with the CFD method based on an efficient dynamic mesh strategy. Ou et al., used the CFD method to analyze the laminar and turbulence flows of a flapping wing robot, and obtained the aerodynamic efficiency of some flapping wing motions [29]. In spite of the fact that some studies have been carried out, the results obtained with the CFD method are not systematic enough, and the results do not include lift and drag coefficient curves, pressure contours and velocity contours at the same time.
Some prototype experiments of flapping wing robots have been performed in recent decades [30,31]. In [32], the aerodynamic forces of a flapping wing robot were measured in experiments, the results showed that the strong thrust can be obtained by increasing the frequency and reducing the flapping amplitude. Harmon et al., placed some reflective markers on a flapping wing in experiments, and used a Vicon measurement system to track the markers for analyzing the wing aerodynamics [33]. In [34], kinematic data in experiments were obtained by tracking small reflective markers on two flapping wings, the data were used to establish flapping membrane wing aerodynamics. In [35], the wing shape of a flapping wing robot was optimized by taking the overall power efficiency as a key objective, and the experiments of a 17.2 g robot were used to verify the conclusion effectiveness. To date, there are still some challenges in the flight experiments. For example, it is difficult to measure all the air flow fields in the robot flights due to the lack of a huge equipment. e experiments are usually carried out outdoors, so the environmental parameters usually have some interference to the experiment results. e complex flight mechanism of a bird-like flapping wing robot is vital for improving its flight efficiency and service life as well as stability. It mainly includes two steps to obtain the complex flight mechanism. e first step is to establish the physical model and the kinematic model as well as the dynamic model of a flapping wing robot. e second step is to conduct simulation experiments and prototype experiments at a low Reynolds number. e structures of this study are as follows: both the physical model and the kinematic model of a bird-like flapping wing robot are established in Section 2. e dynamic model of the robot is established in Section 3. In Section 4, the robot flight is simulated with the CFD method, followed by the prototype experiments. Section 5 and Section 6 give the discussion and conclusions of this study respectively.

e Physical
Model. e physical model of the robot is shown in Figure 2(a). e robot prototype in Figure 2(b) was built with some outsourced components. e prototype includes a body frame, two flapping wings and a tail as well as driving mechanism and transmission mechanism, etc. e aerofoil of the robot is a flat plate. e mass of the robot is nearly 600 g, and the maximum load is 220 g. e wingspan is 1.26 m and the length from head to tail is 0.5 m.
As shown in Figure 3, both the wings and the tail are driven independently. e transmission mechanism of the wings is: motor -gear on the motor-gear A1-gear A2-gear A3-connecting rod-wing bar. Both the gears and the wing bar fix on the body frame with rotation. e above transmission mechanism drives the wings to flap repeatedly with the high accuracy. e wing bar is made of carbon fiber spars, which are with high strength for resisting strong impact forces. e transmission mechanism that drives the tail to swing upwards and downwards is: motor actuator B -connecting rod B1-connecting rod B2-connecting rod B3. e tail swings leftwards and rightwards by controlling motor actuator C. In addition, an electronic speed controller receives signals from remote controller to control the robot flight. Table 1 is the components of the flapping wing robot. Numbers 1.1-1.9 are components related to the wing movements. Numbers 2.1-2.9 are components related to the tail movements. Number 3 is body frame. Numbers 4.1-4.4 are battery or control components.

e Kinematic Model.
e movement of the bird-like flapping wing robot is shown in Figure 4. e wings symmetrically flap with respect to a middle horizontal plane. e wings first flap upwards from position 1 (the horizontal plane) through position 2 to position 3. e wings flap downwards from position 3 through position 4,5,6 to position 7. In the following the wings flap upwards from position 7 through position 8 to position 9 (the horizontal plane).
In this study, the movement of the rigid wings is simplified into only up-down flapping. e kinematic model of the robot is given as follows [19,25]: where α(t) is the flapping angle between the wing and the horizontal plane, A is the flapping amplitude, f is the flapping frequency and t is the time.

The Dynamic Model of the Robot
e Reynolds number of a bird-like flapping wing robot is usually from 5 × 10 3 to 5 × 10 5 during its flights [36]. e RANS equations have been used to simulate complex turbulence of the robot. e average flow can be calculated by RANS equations due to less amount of calculation.
e Spalart-Allmaras model can be used to close the eddy viscosity models in RANS equations. In this study, RANS equations can be written in the Cartesian tensor forms [37,38]: where (a), (b), (c) and (d) are the inertia force, the pressure, the viscous force and the external force acting on the air flow respectively. u is the fluid velocity. ρ denotes the fluid density. p is the pressure. μ denotes the dynamic viscosity. δ is a delta function. −ρu i ′ u j ′ is the Reynolds stresses. e velocity component u i is given as follows: where u i denotes the mean velocity component, and u i ′ is the fluctuating velocity component.
In this study, the Spalart-Allmaras model is used to close the RANS equations at a low Reynolds number. e Spalart-Allmaras model can be expressed as equation (5) [37]: where v denotes the modified turbulent viscosity. P v is the turbulent viscosity production. D v is the turbulent viscosity destruction. Both σ v and C b2 are constants. T v is a userdefined source term.  Journal of Robotics e turbulent viscosity production P v is given as follows [37]: where both C b1 and k are constants. l is the distance from the wall. T denotes a scalar measure of the deformation tensor, and λ ≡ v/v. e mean rate-of-rotation tensor ω ij is modeled as equation (10): e viscous damping function f v1 can be written as equation (11):  where C v1 is a constant. e turbulent viscosity destruction D v is defined as equation (12): where: where C w1 , C w2 and C w3 are constants. C w1 can be calculated by equation (16): e parameter values C b1 , C b2 , C w1 , C w2 , C w3 , C v1 , σ v and k are given in Table 2 [37].
Equations (2)-(4) are RANS equations which have been used to simulate complex turbulence in the robot flights. equations (5) is the Spalart-Allmaras model which is used to close the RANS equations. equation (6)- (16) are the details of some variables in equation (5).
e Reynolds number is defined by equation (17) [12,39]: where L is the chord length of the robot. us, the dynamic model of the robot was established based on the RANS equations and the Spalart-Allmaras model.

The Flight Experiments
Simulation experiments of the robot were performed at a low Reynolds number in Section 4.1. Prototype experiments of the robot were carried out in Section 4.2.

Simulation Experiments.
e steps of the simulation experiments are given as follows: (1) A solid model of the flapping wing robot was built in SOLIDWORKS software, and the file format was saved as.xt or.step. were divided in ANSYS MESHING software on WORKBENCH platform. e moving meshes are required to imitate the movements of the wing. As shown in Figure 5(a), unstructured triangular meshes were established which can simulate the robot movements within a short time, and the file format was saved as.mesh. Details of "mesh" are shown in Figure 5(b) and 5(c).
(4) UDF (User-defined function) codes were programmed based on the kinematic model in Section2. e parameter values in equation (1) were set as follows: the flapping amplitude A � π/6, the flapping frequency f � 50Hz. e fluid velocity u � 5m/s. e Reynolds number was calculated as Re � 7.22 × 10 4 .
(5) e meshes in step (3) were imported to ANSYS FLUENT software, and some options for the robot model were set. In the simulation experiments, it is important to select the viscous model of the robot and set its parameters. As shown in Figure 6(a), the Spalart-Allmaras model was selected, and the parameters in Table 2 were set. Boundary conditions of velocity inlet were set as Figure 6(b). Velocity magnitude was set as 5 m/s, X-component of flow direction was set as 0.9848, Y-component of flow direction was set as 0.17365, specification method of turbulence was selected as modified turbulent viscosity, modified turbulent viscosity was set as 0.0001 m 2 /s. Some other options for the model were also set as follows: the type of solver was set as pressurebased, the velocity formulation of solver was set as absolute, wall motion of both head and tail was selected as stationary wall, shear condition of both head and tail was selected as no slip, roughness constant of both head and tail was set as 0.5. e UDF codes in step (4) were imported in ANSYS FLUENT software.
(6) e position of the flapping wing at each time step was solved for obtaining the results in the transient state and steady state. e post-processing to the results was carried out in ANSYS FLUENT software and TECPLOT software.
Some results of the robot at a low Reynolds number are shown as Figures 7-11.
As shown in Figure 7, the lift and drag coefficient curves are closed as the wing flaps periodically. Most coefficient values are bigger than zero. e integral of the lift or the drag coefficients with respect to the flapping angle is much bigger than zero within a period, which makes the robot fly stably. e upper parts of two curves describe the phase of flapping downwards, which generate the effective aerodynamic forces. Moreover, the aerodynamic force coefficients become bigger while the wing flaps downwards, and they are opposite while the wing flaps upwards. From Figure 7, it can be obtained that the integral of the lift coefficient is bigger than that of the drag coefficient within a period, which indicates the lift is bigger than the drag. In addition, the rate of aerodynamic force curves reaches the biggest value at the upward-downward or downward-upward transition. An air exhaust mechanism should be designed for improving the robot flight efficiency. e curve magnitude in Figure 7 is the  [40,41], but the curve shapes are different with the above references. It expresses the relationship between the aerodynamic force coefficients and the flapping angle in Figure 7 more clearly than that in references [40,41]. It shows that the pressure contours of two wings are symmetrical with respect to the middle plane from Figure 8(a) and 8(b). Most pressure values of the lower surface are higher than zero while those of the upper surface are lower than zero. e closer to the head, the pressure value becomes the lower in Figure 8(b). In Figure 8(c), there are several ellipsoidal pressure vortexes above the upper surface and below the lower surface, because the volume near the wing changed rapidly but the air did not escape in time. ere are the high pressure values at the wing edges that form a small pressure vortex, which is wrapped by several big vortexes. In Figure 8(d), there are two pressure vortex groups with low values below the head connection and tail connection. In Figure 8(e), the pressure value above the wing is lower than that below the wing, which results in an upward lift. As the pressure values are high, the strength of the wing edges and the head connection as well as the tail connection should be increased in order to resist the high Journal of Robotics pressure, which has not been mentioned in previous research results. ere are some similar vorticity contours in reference [24], but it has not mentioned vorticity contours of the robot surfaces, and the specific values of the vorticity contours have not been given in reference [24]. e shapes of pressure iso-surfaces in ere are several arcs in Figure 10(a), and the air velocity below the wing is much lower than that above the wing. It also shows that the air velocity near the upper surface of the wing is low, which is attributed to the air adhesion to the wing. e air velocity gradually becomes lower from the inner layer to the outer layer. e robot is surrounded by velocity contours in Figure 10       back of the wing. e other groups with high values are nearly symmetrical with the above groups with respect to the wing. e centers of the above vortex groups are closed to the wing. ere are also four vortex groups in Figure 10(e) when the wing flaps downwards, but the positions of the vortex groups are opposite to those in Figure 10(d). e robot skin where appears four velocity vortex groups should avoid stitching for preventing air leakages, which has not been proposed in other references. Some velocity contours were given in references [16,42], but there were few figures or lack of specific velocity values. In Figure 11(a), there is a long cylindrical iso-surface above the front of the wing. e robot is wrapped by a sphere iso-surface in Figure 11(b). e robot is covered by a flat iso-surface in Figure 11(c). ere is a flat iso-surface at the back of the wing in Figure 11(d). Figures 11(e)-11(h) are velocity iso-surfaces when the wing flaps upwards. e shape of the iso-surface in Figure 11(e) likes two silkworms, which locate above the front of the robot and below the back of the robot. ere is a very large isosurface in Figure 11

Prototype Experiments.
e robot prototype experiments mainly include the following steps: (1) Before the flight experiments of the robot, the mass center of the robot was adjusted for ensuring the balance between left part and right part, and the tail was kept upwards 30°. (2) When the robot took off, the wing flapped with the highest frequency, and the robot flied against the wind to generate the maximum aerodynamic force for taking off successfully. e Reynolds number is mainly determined by ρ u, L, μ from equation (17). Both the fluid density ρ and the dynamic viscosityμ depend on the external environment. A sunny environment was chosen in the prototype experiments. An anemometer was used to measure wind speed and temperature. e flight velocity of the robotu was kept at 5 m/s by adjusting the remote controller. e chord length of the robot L is a constant. e parameters ρ u, μ had no big deviation in spite of the fact that there were some small changes of the external environment, therefore the robot flied at a low Reynolds number. e remote controller controlled the robot flight by sending signals to the electronic speed controller at the ground station. e motor drove the movements of the wings. e motor was controlled by a handle on the remote controller to realize the robot flight in a straight line. e faster the motor, the faster the robot flied. e motor actuators drove the movements of the tail. e motor actuators were remotely controlled by another handle to realize the robot turning. e faster the motor actuators, the faster the robot turned a corner. e robot flight effects were mainly qualitatively evaluated due to the limitation of experimental conditions. It was considered that the robot flight effects were good when the robot realized stable flight, otherwise it was considered that the flight effects were poor.

Discussion
It is important to obtain the complex flight mechanism at a low Reynolds number for realizing the stable flights of the robot. Some interesting results have been obtained from the flight experiments. For example, most values of aerodynamic force coefficient are bigger than zero that generate the upward lift and forward thrust. e effective aerodynamic forces mainly generate in the phase that the wing flaps downwards instead of flapping upwards. e pressure contours show the big vortexes wrap the small vortexes, which are usually with high pressure values and closed to the wing edges. ere are four velocity vortex groups in total at the front and back of the wing in the velocity contours. With some results in this study, some methods have been obtained for improving the robot flight performances.
Referring to references [19,[37][38][39], the physical model, the kinematic model and the dynamic model were separately established in this study. e 3D physical model with the same parameters as the robot prototype was built in this study, which is much nearer to real robots compared with 1D or 2D model in references [14,39,43,44] e same as references [17,19,40,41], the kinematic model of the robot was established in this study. e dynamic model was established based on the RANS equations and the Spalart-Allmaras model in this study, which can effectively describe the air flow fields around the robot. e strip method was used to establish the dynamic model of flapping wing robots in references [19,34,40,41,45,46], by which the aerodynamic forces have been approximately calculated, but the results are less accurate than that in this study. In references [20,43,47,48], the dynamic model of flapping wing robot was established with the moment equilibrium method, which can approximately analyze the complex dynamic characteristics, but it is difficult to simulate the air flow fields around flapping wings in realities. In summary, both the physical model and the dynamic model are developed in this study compared with the previous research results.
e simulation experiments were carried out in this study. e lift and drag coefficient curves within a period were obtained to describe the aerodynamic forces, both the pressure contours and the velocity contours were used to analyze the complex flight mechanism in this study. Some simulation results in this study are similar to those in reference [24]. e aerodynamic forces of the robot were obtained with some non fluid simulation methods in references [19,34,40,41,45,46,49], which have not described the air flow fields and been less accurate than the results in this study. In addition, some experiments were performed in references [16,42,50] but which did not study steady situations and transient situations as well as prototype experiments together.
e simulation experiments were explored more in this study compared with the previous research results. e prototype experiments were carried out in this study similar to references [32,33,35], and the robot flied stably at a low Reynolds number. However, the air flow fields around the robot have not been measured in the prototype experiments, and the robot has not flied in some complex situations, e.g., severe weathers. Some methods for improving the robot flight performances have been proposed in this study, which have not been given in the previous research results.
In spite of the fact that the complex flight mechanism of the robot has been analyzed in this study, the following challenges should be further explored: (1) A bird-like robot prototype should be further optimized with bionics and biomimetics, and a more complex kinematic model includes the relative incidence angle should be established for imitating the movement of real birds. (2) Some flight experiments should be improved with advanced instruments and methods, e.g., a camera measurement or a wind tunnel should be used to measure the air flow fields in prototype experiments. A quantitative measurement method should be adopted to accurately evaluate the robot flight effects. (3) e complex flight mechanism of the robot under severe weathers should be further studied in the future.

Conclusions
e complex flight mechanism of a bird-like flapping wing robot at a low Reynolds number was studied in this study. Both the physical model and the kinematic model of the robot were first built, and the dynamic model was also established based on the RANS equations and the Spalart-Allmaras turbulence model. e flight experiments of the robot were performed. Some important results have been obtained from the flight experiments: (1) e robot flies upwards and forwards mainly because aerodynamic force coefficient is bigger than zero within a period. It is the upward lift and forward thrust when the wings flap downwards. e rate of the coefficient curves reaches the biggest value when the flapping direction changes. (2) ere are several pressure vortexes near to the flapping wing, the big vortexes wrapped the small vortexes, which are with the high pressure values and usually appear at the wing edges. e above pressure vortexes make the robot fly forwards. (3) It indicates that the air velocity below the wing is much lower than that above the wing in the steady state. When the wing flaps downwards or upwards, there are four velocity vortex groups in the transient state.
Some methods for improving the robot flight performances have been obtained from the flight experiments. (1) As the upward or downward lift mainly generates in the phase that the wing flaps downwards or upwards, an air exhaust mechanism should be designed for improving the robot flight efficiency. (2) e strength of the wing edges and the head connection as well as the tail connection should be increased for resisting the high pressure. (3) e robot skin where appears four velocity vortex groups should avoid stitching in order to prevent air leakages.
Data Availability e figures and estimation results data used to support the findings of this study are included within the article.

Conflicts of Interest
e authors declare that there is no conflict of interest regarding the publication of this paper.