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

Authors: Changtao Ding; Xiating Yao; Chengyao Liu
Journal: Journal of Robotics
 (2022)

---

## Abstract

The flight mechanism of a bird-like flapping 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 first 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 flight experiments were carried out and the results were discussed. Lift and drag coefficient curves show that it generates upward lift and forward thrust in the phase that the wing flaps downwards, the rate of the coefficient curves is the biggest when the flapping 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 flight efficiency 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 flights of the flapping wing robot with the high efficiency.


---

## Body

## 1. 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–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–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].Figure 1
Several flapping wing robots. where (a) :the robot designed by Chirarattananon; (b) :the robot designed by Leys; (c) :the robot designed by Grand.
(a)(b)(c)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. The 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. Therefore, it is necessary to establish a dynamic model that is closer to the real flight for obtaining the robot flight mechanism.The 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. The experiments are usually carried out outdoors, so the environmental parameters usually have some interference to the experiment results.The 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. The first step is to establish the physical model and the kinematic model as well as the dynamic model of a flapping wing robot. The second step is to conduct simulation experiments and prototype experiments at a low Reynolds number. The 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 Section2. The 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.


## 2. The Physical Model and the Kinematic Model

### 2.1. The Physical Model

The physical model of the robot is shown in Figure2(a). The robot prototype in Figure 2(b) was built with some outsourced components. The prototype includes a body frame, two flapping wings and a tail as well as driving mechanism and transmission mechanism, etc. The aerofoil of the robot is a flat plate. The mass of the robot is nearly 600 g, and the maximum load is 220 g. The wingspan is 1.26 m and the length from head to tail is 0.5 m.Figure 2
A bird-like flapping wing robot. where (a): the physical model of the robot; (b): the robot prototype.
(a)(b)As shown in Figure3, both the wings and the tail are driven independently. The 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. The above transmission mechanism drives the wings to flap repeatedly with the high accuracy. The wing bar is made of carbon fiber spars, which are with high strength for resisting strong impact forces. The transmission mechanism that drives the tail to swing upwards and downwards is: motor actuator B - connecting rod B1- connecting rod B2- connecting rod B3. The 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.Figure 3
The components of the flapping wing robot. where (a): the robot in global view; (b): the robot in bottom view.
(a)(b)Table1 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.Table 1
The components of the flapping wing robot.
NumberComponentQuantityDescriptions1.1Motor1It drives the movement of two flapping wings.1.2Gear on the motor1It coaxially fixes with motor and meshes with gear A1.1.3Gear A11It meshes with gear on the motor and coaxially fixes with gear A2.1.4Gear A22It coaxially fixes with gear A1 and meshes with gear A3.1.5Gear A32It meshes with gear A2 and rotationally connects with connecting rod.1.6Connecting rod2It rotationally connects with gear A3 at one end, and it rotationally connects with wing bar by a ball joint at the other end.1.7Wing bar2It rotationally connects with connecting rod and fixes with ribs in flapping wing.1.8Rib in flapping wing10It fixes at wing bar to prop up skin on flapping wing.1.9Skin on flapping wing1It is propped up by wing bars and ribs in flapping wing.2.1Motor actuator B1It fixes in body frame and drives the tail to swing upwards and downwards.2.2Connecting rod B11It rotationally connects with motor actuator B and connecting rod B2.2.3Connecting rod B21It rotationally connects with connecting rod B1 and connecting rod B3.2.4Connecting rod B31It rotationally connects with motor actuator B2 and connecting frame.2.5Connecting frame1It rotationally connects with connecting rod B3 and body frame.2.6Motor actuator C1It fixes in connecting frame and drives the tail to swing leftwards and rightwards.2.7Tail plate1It rotationally connects with motor actuator C, and fixes ribs in tail.2.8Rib in tail4It fixes at tail plate to prop up skin on tail.2.9Skin on tail1It is propped up by tail plate and ribs in tail.3Body frame1It is the robot body that fixes components.4.1Battery1It provides power for motor and motor actuator B as well as motor actuator C.4.2Electronic speed controller1It connects with battery and 12 bits 6 channel receiver as well as motor.4.312 bits 6 channel receiver1It connects with electronic speed controller and motor actuator B as well as motor actuator C.4.4Remote controller1It remotely sends signals to control the robot movements.


### 2.2. The Kinematic Model

The movement of the bird-like flapping wing robot is shown in Figure4. The wings symmetrically flap with respect to a middle horizontal plane. The wings first flap upwards from position 1 (the horizontal plane) through position 2 to position 3. The 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).Figure 4
The movement of the bird-like flapping wing robot.In this study, the movement of the rigid wings is simplified into only up-down flapping. The kinematic model of the robot is given as follows [19, 25]:(1)αt=Asin2πft,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.


## 2.1. The Physical Model

The physical model of the robot is shown in Figure2(a). The robot prototype in Figure 2(b) was built with some outsourced components. The prototype includes a body frame, two flapping wings and a tail as well as driving mechanism and transmission mechanism, etc. The aerofoil of the robot is a flat plate. The mass of the robot is nearly 600 g, and the maximum load is 220 g. The wingspan is 1.26 m and the length from head to tail is 0.5 m.Figure 2
A bird-like flapping wing robot. where (a): the physical model of the robot; (b): the robot prototype.
(a)(b)As shown in Figure3, both the wings and the tail are driven independently. The 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. The above transmission mechanism drives the wings to flap repeatedly with the high accuracy. The wing bar is made of carbon fiber spars, which are with high strength for resisting strong impact forces. The transmission mechanism that drives the tail to swing upwards and downwards is: motor actuator B - connecting rod B1- connecting rod B2- connecting rod B3. The 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.Figure 3
The components of the flapping wing robot. where (a): the robot in global view; (b): the robot in bottom view.
(a)(b)Table1 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.Table 1
The components of the flapping wing robot.
NumberComponentQuantityDescriptions1.1Motor1It drives the movement of two flapping wings.1.2Gear on the motor1It coaxially fixes with motor and meshes with gear A1.1.3Gear A11It meshes with gear on the motor and coaxially fixes with gear A2.1.4Gear A22It coaxially fixes with gear A1 and meshes with gear A3.1.5Gear A32It meshes with gear A2 and rotationally connects with connecting rod.1.6Connecting rod2It rotationally connects with gear A3 at one end, and it rotationally connects with wing bar by a ball joint at the other end.1.7Wing bar2It rotationally connects with connecting rod and fixes with ribs in flapping wing.1.8Rib in flapping wing10It fixes at wing bar to prop up skin on flapping wing.1.9Skin on flapping wing1It is propped up by wing bars and ribs in flapping wing.2.1Motor actuator B1It fixes in body frame and drives the tail to swing upwards and downwards.2.2Connecting rod B11It rotationally connects with motor actuator B and connecting rod B2.2.3Connecting rod B21It rotationally connects with connecting rod B1 and connecting rod B3.2.4Connecting rod B31It rotationally connects with motor actuator B2 and connecting frame.2.5Connecting frame1It rotationally connects with connecting rod B3 and body frame.2.6Motor actuator C1It fixes in connecting frame and drives the tail to swing leftwards and rightwards.2.7Tail plate1It rotationally connects with motor actuator C, and fixes ribs in tail.2.8Rib in tail4It fixes at tail plate to prop up skin on tail.2.9Skin on tail1It is propped up by tail plate and ribs in tail.3Body frame1It is the robot body that fixes components.4.1Battery1It provides power for motor and motor actuator B as well as motor actuator C.4.2Electronic speed controller1It connects with battery and 12 bits 6 channel receiver as well as motor.4.312 bits 6 channel receiver1It connects with electronic speed controller and motor actuator B as well as motor actuator C.4.4Remote controller1It remotely sends signals to control the robot movements.


## 2.2. The Kinematic Model

The movement of the bird-like flapping wing robot is shown in Figure4. The wings symmetrically flap with respect to a middle horizontal plane. The wings first flap upwards from position 1 (the horizontal plane) through position 2 to position 3. The 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).Figure 4
The movement of the bird-like flapping wing robot.In this study, the movement of the rigid wings is simplified into only up-down flapping. The kinematic model of the robot is given as follows [19, 25]:(1)αt=Asin2πft,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.


## 3. The Dynamic Model of the Robot

The Reynolds number of a bird-like flapping wing robot is usually from5×103 to 5×105 during its flights [36]. The RANS equations have been used to simulate complex turbulence of the robot. The average flow can be calculated by RANS equations due to less amount of calculation. The 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]:(2)∂ρ∂t+∂∂xiρui=0,(3)∂∂tρui+∂∂xjρuiuj︸a=−∂p∂xi︸b+∂∂xjμ∂ui∂xj+∂uj∂xi−23δij∂u1∂x1︸c+∂∂xj−ρui′uj′¯︸d,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. −ρui′uj′¯ is the Reynolds stresses.The velocity componentui is given as follows:(4)ui=u¯i+ui′i=1,2,3,where u¯i denotes the mean velocity component, and ui′ is the fluctuating velocity component.In this study, the Spalart-Allmaras model is used to close the RANS equations at a low Reynolds number. The Spalart-Allmaras model can be expressed as equation (5) [37]:(5)∂∂tρv˜+∂∂xiρv˜ui=Pv+1σv˜∂∂xjμ+ρv˜∂v˜∂xj+Cb2ρ∂v˜∂xj2−Dv+Tv˜,where v˜ denotes the modified turbulent viscosity. Pv is the turbulent viscosity production. Dv is the turbulent viscosity destruction. Both σv˜ and Cb2 are constants. Tv˜ is a user-defined source term.The turbulent viscosity productionPv is given as follows [37]:(6)Pv=Cb1ρT˜v˜,(7)T˜≡T+v˜k2l2fv2,(8)T≡2ωijωij,(9)fv2≡1−λ1+λfv1,where both Cb1 and k are constants. l is the distance from the wall. T denotes a scalar measure of the deformation tensor, and λ≡v˜/v.The mean rate-of-rotation tensorωij is modeled as equation (10):(10)ωij=12∂ui∂xj−∂uj∂xi.The viscous damping functionfv1 can be written as equation (11):(11)fv1=λ3λ3+Cv13,where Cv1 is a constant.The turbulent viscosity destructionDv is defined as equation (12):(12)Dv=Cw1ρfwv˜l2,where:(13)fw=g1+Cw36g6+Cw361/6,(14)g=r+Cw2r6−r,(15)r≡v˜T˜k2l2,where Cw1, Cw2 and Cw3 are constants.Cw1 can be calculated by equation (16):(16)Cw1=Cb1k2+1+Cb2σv˜.The parameter valuesCb1, Cb2, Cw1, Cw2, Cw3, Cv1, σv˜ and k are given in Table 2 [37].Table 2
Some parameter values.
ParametersCb1Cb2Cw1Cw2Cw3Cv1σv˜kValues0.13550.6223.2060.32.07.12/30.4187Equations (2)–(4) are RANS equations which have been used to simulate complex turbulence in the robot flights. Equation (5) is the Spalart-Allmaras model which is used to close the RANS equations. Equations (6)–(16) are the details of some variables in equation (5).The Reynolds number is defined by equation (17) [12, 39]:(17)Re=ρuLμ,where L is the chord length of the robot.Thus, the dynamic model of the robot was established based on the RANS equations and the Spalart-Allmaras model.


## 4. The Flight Experiments

Simulation experiments of the robot were performed at a low Reynolds number in Section4.1. Prototype experiments of the robot were carried out in Section 4.2.


### 4.1. Simulation Experiments

The 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.(2)
The model file in step (1) was imported to ANSYS DM software on WORKBENCH platform, and the calculation domains of both the robot and the air flow fields were set based on real flight situations.(3)
The meshes of both the robot and the air flow fields were divided in ANSYS MESHING software on WORKBENCH platform. The moving meshes are required to imitate the movements of the wing. As shown in Figure5(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).Figure 5
The meshes of both the robot and the air flow fields. (a): unstructured triangular meshes; (b) and (c): details of “mesh”.
(a)(b)(c)(4) UDF (User-defined function) codes were programmed based on the kinematic model in Section2. The parameter values in equation (1) were set as follows: the flapping amplitude A=π/6, the flapping frequency f=50Hz. The fluid velocity u=5m/s. The Reynolds number was calculated as Re=7.22×104.(5) The 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 Figure6(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 m2/s. Some other options for the model were also set as follows: the type of solver was set as pressure-based, 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. The UDF codes in step (4) were imported in ANSYS FLUENT software.Figure 6
The settings of viscous model and velocity inlet in ANSYS FLUENT software. (a): the settings of viscous model; (b): the settings of velocity inlet.
(a)(b)(6) The position of the flapping wing at each time step was solved for obtaining the results in the transient state and steady state. The 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 Figures7–11.Figure 7
Lift and drag coefficient curves of the bird-like flapping wing robot.Figure 8
Pressure contours of the flapping wing robot. (a) and (b): the lower and upper surfaces of the robot in the steady state; (c): 3D view along the symmetrical plane in the steady state; (d): 3D view in the transient state when the wing flaps upwards; (e): 3D view in the transient state when the wing flaps downwards.
(a)(b)(c)(d)(e)Figure 9
Pressure iso-surfaces of the flapping wing robot. (a), (b), (c) and (d):static pressureP = −4, −1, 1, and 4 (Pa) in the steady static; (e), (f), (g) and (h): static pressure P = 0, 1200, 1850 and 2000 (Pa) when the wing flaps upwards; (i), (j), (k) and (l): static pressure P = 0, 1000, 2000 and 3000 (Pa) when the wing flaps downwards.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)Figure 10
Velocity contours of the flapping wing robot. (a), (b) and (c): right view, top view and front view of the robot in the steady state; (d): front view in the transient state when the wing flaps upwards; (e): front view in the transient state when the wing flaps downwards.
(a)(b)(c)(d)(e)Figure 11
Velocity iso-surfaces of the flapping wing robot. (a), (b), (c) and (d): air velocityv = 6, 5, 4.8 and 3 in the steady state; (e), (f), (g) and (h): air velocity v = 0, 5, 6 and 10 in the transient state when the wing flaps upwards; (i), (j), (k) and (l): air velocity v = 0, 4, 6 and 10 in the transient state when the wing flaps downwards.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)As shown in Figure7, the lift and drag coefficient curves are closed as the wing flaps periodically. Most coefficient values are bigger than zero. The 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. The 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. The curve magnitude in Figure 7 is the same as references [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 Figure8(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. The 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. There 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 pressure, which has not been mentioned in previous research results. There 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].The shapes of pressure iso-surfaces in Figures9 are bilaterally symmetrical with respect to the middle plane. The pressure iso-surfaces in Figures 9(a) and 9(b) are similar to the smooth cylinders that are above the flapping wing. The iso-surfaces in Figures 9(a) and 9(b) are nearly symmetrical with those in Figures 9(d) and 9(c). The pressure iso-surfaces in Figures 9(e)–9(h) are nearly symmetrical with those in Figures 9(i)–9(l) due to the opposite flapping direction. The pressure iso-surfaces in Figures 9(e) and 9(i) resemble the ring shapes, and the iso-surfaces in Figures 9(f)–9(h) and 9(j)–9(l) resemble the spherical shapes.There are several arcs in Figure10(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. The air velocity gradually becomes lower from the inner layer to the outer layer. The robot is surrounded by velocity contours in Figure 10(b). From Figure 10(c), three velocity vortexes with high values appear above the front of the wing, two velocity vortexes with low values locate below the front of the wing. There are four velocity vortex groups in Figure 10(d) when the wing flaps upwards. Two groups with low values locate above the front of the wing and below the back of the wing. The other groups with high values are nearly symmetrical with the above groups with respect to the wing. The centers of the above vortex groups are closed to the wing. There 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). The 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 Figure11(a), there is a long cylindrical iso-surface above the front of the wing. The robot is wrapped by a sphere iso-surface in Figure 11(b). The robot is covered by a flat iso-surface in Figure 11(c). There 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. The 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. There is a very large iso-surface in Figure 11(f). The iso-surfaces in Figures 11(g) and 11(h) also like silkworms, which locate above the back of the robot and below the front of the robot. The velocity iso-surfaces in Figures 11(i)–11(l) like silkworms, which appear when the wing flaps downwards. The iso-surfaces in Figures 11(i) and 11(j) locate above the back of the robot and below the front of the robot. The iso-surfaces in Figures 11(k) and 11(l) locate above the front of the robot and below the back of the robot.


### 4.2. Prototype Experiments

The 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. (3) When the robot flied stably in the sky, the flapping frequency of the wing turned to 1/2 or 2/3 of that in taking off phase. The robot flied towards different direction by adjusting its tail. The robot would turn left or right when the tail turned left or right. (4) The robot floated down freely at the end of the flight for avoiding severe impacts with the ground.The flight experiments of the robot are shown in Figure12, the robot flied stably and freely in the sky at a low Reynolds number. During the flight, the robot repeatedly flapped its wings to generate enough aerodynamic forces, and the robot turned left or right by turning its tail left or right.Figure 12
Prototype experiments of the robot.The 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. The flight velocity of the robotu was kept at 5 m/s by adjusting the remote controller. The chord length of the robot L is a constant. The 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.The remote controller controlled the robot flight by sending signals to the electronic speed controller at the ground station. The motor drove the movements of the wings. The motor was controlled by a handle on the remote controller to realize the robot flight in a straight line. The faster the motor, the faster the robot flied. The motor actuators drove the movements of the tail. The motor actuators were remotely controlled by another handle to realize the robot turning. The faster the motor actuators, the faster the robot turned a corner.The 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.


## 4.1. Simulation Experiments

The 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.(2)
The model file in step (1) was imported to ANSYS DM software on WORKBENCH platform, and the calculation domains of both the robot and the air flow fields were set based on real flight situations.(3)
The meshes of both the robot and the air flow fields were divided in ANSYS MESHING software on WORKBENCH platform. The moving meshes are required to imitate the movements of the wing. As shown in Figure5(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).Figure 5
The meshes of both the robot and the air flow fields. (a): unstructured triangular meshes; (b) and (c): details of “mesh”.
(a)(b)(c)(4) UDF (User-defined function) codes were programmed based on the kinematic model in Section2. The parameter values in equation (1) were set as follows: the flapping amplitude A=π/6, the flapping frequency f=50Hz. The fluid velocity u=5m/s. The Reynolds number was calculated as Re=7.22×104.(5) The 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 Figure6(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 m2/s. Some other options for the model were also set as follows: the type of solver was set as pressure-based, 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. The UDF codes in step (4) were imported in ANSYS FLUENT software.Figure 6
The settings of viscous model and velocity inlet in ANSYS FLUENT software. (a): the settings of viscous model; (b): the settings of velocity inlet.
(a)(b)(6) The position of the flapping wing at each time step was solved for obtaining the results in the transient state and steady state. The 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 Figures7–11.Figure 7
Lift and drag coefficient curves of the bird-like flapping wing robot.Figure 8
Pressure contours of the flapping wing robot. (a) and (b): the lower and upper surfaces of the robot in the steady state; (c): 3D view along the symmetrical plane in the steady state; (d): 3D view in the transient state when the wing flaps upwards; (e): 3D view in the transient state when the wing flaps downwards.
(a)(b)(c)(d)(e)Figure 9
Pressure iso-surfaces of the flapping wing robot. (a), (b), (c) and (d):static pressureP = −4, −1, 1, and 4 (Pa) in the steady static; (e), (f), (g) and (h): static pressure P = 0, 1200, 1850 and 2000 (Pa) when the wing flaps upwards; (i), (j), (k) and (l): static pressure P = 0, 1000, 2000 and 3000 (Pa) when the wing flaps downwards.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)Figure 10
Velocity contours of the flapping wing robot. (a), (b) and (c): right view, top view and front view of the robot in the steady state; (d): front view in the transient state when the wing flaps upwards; (e): front view in the transient state when the wing flaps downwards.
(a)(b)(c)(d)(e)Figure 11
Velocity iso-surfaces of the flapping wing robot. (a), (b), (c) and (d): air velocityv = 6, 5, 4.8 and 3 in the steady state; (e), (f), (g) and (h): air velocity v = 0, 5, 6 and 10 in the transient state when the wing flaps upwards; (i), (j), (k) and (l): air velocity v = 0, 4, 6 and 10 in the transient state when the wing flaps downwards.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)As shown in Figure7, the lift and drag coefficient curves are closed as the wing flaps periodically. Most coefficient values are bigger than zero. The 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. The 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. The curve magnitude in Figure 7 is the same as references [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 Figure8(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. The 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. There 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 pressure, which has not been mentioned in previous research results. There 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].The shapes of pressure iso-surfaces in Figures9 are bilaterally symmetrical with respect to the middle plane. The pressure iso-surfaces in Figures 9(a) and 9(b) are similar to the smooth cylinders that are above the flapping wing. The iso-surfaces in Figures 9(a) and 9(b) are nearly symmetrical with those in Figures 9(d) and 9(c). The pressure iso-surfaces in Figures 9(e)–9(h) are nearly symmetrical with those in Figures 9(i)–9(l) due to the opposite flapping direction. The pressure iso-surfaces in Figures 9(e) and 9(i) resemble the ring shapes, and the iso-surfaces in Figures 9(f)–9(h) and 9(j)–9(l) resemble the spherical shapes.There are several arcs in Figure10(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. The air velocity gradually becomes lower from the inner layer to the outer layer. The robot is surrounded by velocity contours in Figure 10(b). From Figure 10(c), three velocity vortexes with high values appear above the front of the wing, two velocity vortexes with low values locate below the front of the wing. There are four velocity vortex groups in Figure 10(d) when the wing flaps upwards. Two groups with low values locate above the front of the wing and below the back of the wing. The other groups with high values are nearly symmetrical with the above groups with respect to the wing. The centers of the above vortex groups are closed to the wing. There 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). The 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 Figure11(a), there is a long cylindrical iso-surface above the front of the wing. The robot is wrapped by a sphere iso-surface in Figure 11(b). The robot is covered by a flat iso-surface in Figure 11(c). There 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. The 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. There is a very large iso-surface in Figure 11(f). The iso-surfaces in Figures 11(g) and 11(h) also like silkworms, which locate above the back of the robot and below the front of the robot. The velocity iso-surfaces in Figures 11(i)–11(l) like silkworms, which appear when the wing flaps downwards. The iso-surfaces in Figures 11(i) and 11(j) locate above the back of the robot and below the front of the robot. The iso-surfaces in Figures 11(k) and 11(l) locate above the front of the robot and below the back of the robot.


## 4.2. Prototype Experiments

The 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. (3) When the robot flied stably in the sky, the flapping frequency of the wing turned to 1/2 or 2/3 of that in taking off phase. The robot flied towards different direction by adjusting its tail. The robot would turn left or right when the tail turned left or right. (4) The robot floated down freely at the end of the flight for avoiding severe impacts with the ground.The flight experiments of the robot are shown in Figure12, the robot flied stably and freely in the sky at a low Reynolds number. During the flight, the robot repeatedly flapped its wings to generate enough aerodynamic forces, and the robot turned left or right by turning its tail left or right.Figure 12
Prototype experiments of the robot.The 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. The flight velocity of the robotu was kept at 5 m/s by adjusting the remote controller. The chord length of the robot L is a constant. The 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.The remote controller controlled the robot flight by sending signals to the electronic speed controller at the ground station. The motor drove the movements of the wings. The motor was controlled by a handle on the remote controller to realize the robot flight in a straight line. The faster the motor, the faster the robot flied. The motor actuators drove the movements of the tail. The motor actuators were remotely controlled by another handle to realize the robot turning. The faster the motor actuators, the faster the robot turned a corner.The 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.


## 5. 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. The effective aerodynamic forces mainly generate in the phase that the wing flaps downwards instead of flapping upwards. The pressure contours show the big vortexes wrap the small vortexes, which are usually with high pressure values and closed to the wing edges. There 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–39], the physical model, the kinematic model and the dynamic model were separately established in this study. The 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] The same as references [17, 19, 40, 41], the kinematic model of the robot was established in this study. The 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. The 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.The simulation experiments were carried out in this study. The 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]. The 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. The simulation experiments were explored more in this study compared with the previous research results.The 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)The complex flight mechanism of the robot under severe weathers should be further studied in the future.


## 6. Conclusions

The 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.The flight experiments of the robot were performed. Some important results have been obtained from the flight experiments: (1) The 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. The rate of the coefficient curves reaches the biggest value when the flapping direction changes. (2) There 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. The 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) The strength of the wing edges and the head connection as well as the tail connection should be increased for resisting the high pressure. (3) The robot skin where appears four velocity vortex groups should avoid stitching in order to prevent air leakages.



---

*Source: 6638104-2022-04-18.xml*