A Linearized Model for an Ornithopter in Gliding Flight: Experiments and Simulations

This work studies the accuracy of a simple but effective analytical model for a flapping-wings UAV in longitudinal gliding flight configuration comparing it with experimental results of a real ornithopter. The aerodynamic forces are modeled following the linearized potential theory for a flat plate in gliding configuration, extended to flappingwing episodes modeled also by the (now unsteady) linear potential theory, which are studied numerically. In the gliding configuration, the model reaches a steady-state descent at given terminal velocity and pitching and gliding angles, governed by the wings and tail position. In the flapping-wing configuration, it is noticed that the vehicle can increase its flight velocity and perform climbing episodes. A realistic simulation tool based on Unreal Engine 4 was developed to visualize the effect of the tail position and flapping frequencies and amplitudes on the ornithopter flight in real time. The paper also includes the experimental validation of the gliding flight and the data has been released for the community.


I. INTRODUCTION
The development of bio-inspired UAVs alternating gliding and flapping-wing episodes is a quite challenging topic nowadays. Contrary to multi-rotor vehicles, ornithopters can save a large amount of energy by simply extending their wings and exploiting thermal currents as real birds do in soaring flights. The role of nature in this framework has been decisive because avian flight has been evolved to achieve very low gliding angles. This fact is due to a simple force equilibrium: roughly speaking, the lift force compensates the weight, and the thrust force the drag. The vertical-tohorizontal velocity ratio is of course prescribed for each specie, being it a measurement of the flight performance. For instance, Tucker and Parrot [1] found that the gliding angle of falcons is in the range 5.5 • − 8 • , reporting as well those of pigeons (12 • − 15 • ), gulls (7 • − 11.5 • ) and vultures (3.5 • − 5 • ), and comparing all of them against the performance of a sailplane (1.5 • − 2.47 • ). Alternatively, the gliding efficiency can also be expressed in terms of the Liftto-Drag ratio L/D, which is also a popular scale in avian flight [2].
But gliding birds also need to flap their wings occasionally either to increase their flight velocity or simply to recover the height invested in the previous glide episode (intermittent flight). This flapping flight configuration is characterized by specific frequencies and amplitudes so that most species (see e.g. [3]- [5]) travel in the interval 0.2 < St = f A/U < 0.4 (f is the flapping frequency, A the stroke amplitude and U the forward velocity), according to [6]. In this interval, it is found the maximum propulsive efficiency [7]. On the other hand, the unsteady lift originated by a twodimensional flapping airfoil in steady flight was first modeled mathematically by Theodorsen [8] in 1935, extending Garrick [9] the study for the thrust force one year later. These models are based on the linearized potential theory, providing quite accurate results for sufficiently high Reynolds numbers with low flapping amplitudes at moderate frequencies, i.e. when separation phenomena does not occur. Very recently, the work of Fernandez-Feria [10] reformulated Garrick's work from the perspective of the impulse theory, improving the agreement of thrust force and efficiency [11].
Thus, the main objective of this work is the validation of a theoretical model based on the two-dimensional formulation of [8] and [10], only accounting for gliding flight, and corrected to take into account three-dimensional effects associated to the finite span of actual wings.
Additionally, the realistic simulator Unreal Engine 4 (UE4) [12] has been used to allow the user to visualize instantaneously the behavior of the model by acting upon the tail deflection and exploring the effect of the frequencies and amplitudes in flapping-wing flight episodes.
The paper is organized as follows. First of all, the theoretical model and the meaning of the parameters are described in Section II. Then, in Section III, a set of experiments are shown to validate the proposed model. Section IV discourses about the numerical validation results and it presents a realistic environment simulator based on Unreal Engine using the proposed model in order to provide a framework to develop algorithms for ornithopters. Finally, Section V summarizes the contributions of this article and some future steps are drawn.

II. THEORETICAL MODEL A. Nondimensional Newton-Euler equations
Let {U c , L c , t c } be a vector defining the characteristic velocity, length and time of the problem so that with m the mass, g the gravitational acceleration, c the mean chord length, ρ the air density and A the wing aspect ratio. Thus, the set of nondimensional equations describing the dynamic behavior of a point-mass ornithopter in longitudinal flight are given by where u and w are the velocity components in the body axis, U b is their magnitude, θ the pitch angle and ω its temporal derivative to reduce the order of the system. The coefficients C L and C Di stand for the lift and the induced drag generated by wings (those relative to the tail incorporate the subscript t), and Li indicates the Lighthill number [13], i.e. a scaling of the parasitic drag [14], Additionally, the following nondimensional parameters appear for the relation between the tail and wing areas, the relation between the producing-moment distances from the center of mass to the tail and the wings (see Fig. 1), the dimensionless mass, and the nondimensional inverse of the moment of inertia, respectively. The system (2)-(5) is solved analytically in steady-state by perturbation methods around a typical gliding flight configuration (β π/2 and θ π/2). It yields from [15] with θ i a second order correction accounting for the induced drag effect, Alternatively, one can write the solutions (8)- (11) in polar form as being γ the so-called gliding angle.

B. Modeling of the aerodynamic forces for gliding and flapping-wings flight configurations
The lift coefficients are modeled according to the potential theory, being the wing considered as a flat plate [16] and the tail as a delta flat plate [17], with α and α t the angles of attack of the wing and the tail, respectively. Note that both coefficients are properly corrected to allow for finite span lifting surfaces. From the previous definitions, the induced drag coefficients are derived as In order to avoid unphysical predictions in stall and poststall zones, the trend of (15)-(16) was saturated to α = 15 • and α t = 35 • (see [17] for a more detailed explanation).
On the other hand, birds in gliding flight also flap their wings occasionally to accelerate themselves controlling their flight path. The flapping flight configuration is essentially characterized by the reduced frequency k = 2πf /U b and the nondimensional amplitude h 0 (here c = 2 because the lengths are scaled with c/2). For a wing undergoing a harmonic heaving motion the lift coefficient is defined by Theodorsen [8], [10], as where F (k) and G(k) are the real and the imaginary part of Theodorsen's function C(k) = F (k) + iG(k). Note that a clear separation was made in (19) to consider the finite span of the wing. Following [18], [19], the factor A/(A + 2) is intended for the correction of circulatory terms, and A/(A + 1) for the added mass. Similarly, the thrust force originated by a flapping wing [10], [11], is given by where C T L is the projection of the lift force onto the velocity direction, C T1 is the added-mass contribution (not included in this case because there is no oscillatory pitching motion) and C T2 accounts for the wing-wake interaction (circulatory contribution). These contributions are defined as with F 1 (k) and G 1 (k) the real and the imaginary parts of the function C 1 (k) = F 1 (k) + iG 1 (k) defined in [10].
Finally, note that the formulations in [8], [10], [11] converge to the well-known solutions ∼ 2πα and ∼ 2πα 2 for the respective lift and thrust force of a flat plate when kh 0 → 0. To conclude, it is important to note that the aerodynamic effect of the flapping wings is only explored in the realistic simulation tool of Sec. IV, and the actual experimental validation is against the gliding theoretical model.

III. EXPERIMENTAL SETUP A. Ornithopter device
A commercially available ornithopter frame from the company CarbonSail was chosen. Its fuselage is practically plane and made of carbon fiber, with two connected flapping wings and a triangular flat plate. The wings and the tail are made from nylon ripstop and are extended over the main structure made from carbon fiber rods. For the performance of the flapping movement, there is an actuation mechanism consisting of two parts. The first one is made of machined aluminum and converts the gearbox spinning motion to the alternative vertical performed by the wings, and the second one is a gearbox made up of a set of gears of brass and acetal.
The ornithopter dimensions are indicated in Fig. 2 and Table I, corresponding these dimensional magnitudes to the  nondimensional parameters in Table II.

B. Hardware specifications
In this subsection, it is presented the hardware required to perform the experimental measurements of the ornithopter.
An on board lightweight hardware system has been developed to obtain appropriate flight conditions (note that the velocity components, and of course the gliding angle, depends on the mass through the scaling in (1)). The devices installed in the ornithopter consist of a Raspberry Pi Zero, an Arduino MICRO board, and an Inertial Measurement Unit (POLOLU AltIMU-10 V5). The latter includes an accelerometer, a compass, a gyroscope and an altimeter. The sensitivity range of the accelerometer is ±2g and the range of the gyroscope is ±2000º/s. Both devices have a 204 Hz refresh rate and 16-bit reading per axis. On the other hand, the micro-board controls the ornithopter actuators and save the odometry data from the IMU. The Raspberry is intended to save the data provided from the IMU, and the Arduino selects the control mode. However, the radio tail control allows the pilot to manipulate the tail deflection in flight in order to ensure stability and safety. The system is powered by a 2S battery with 500 mAh. Figure 3 displays the hardware connections between all the devices. It is important to note that in order to overcome the processing delay introduced by the IMU into the Arduino, the Raspberry Pi was placed to acquire and process the measures. Thus, the experimental data are acquired from the on board system and three cameras tracking the platform as shown in Fig. 6. All camera poses are obtained from a Theodolite Total Station (TS) in a local coordinate system with high accuracy. This setup provides a method to compute the trajectory without the drift caused by the IMU.
Furthermore, in order to control the flight initial conditions, it was essential to develop the launcher platform. This guarantees the repetitiveness and post-processing purposes and even determines the transient period until the steady state is achieved. The platform developed is shown in Fig. 4. Structurally, this platform consists of three different parts: the sliding, the propulsion, and the launcher system. The sliding system involves two steel bars of 8 mm and a sliding roller over the bar. The propulsion system is driven by an electric motor connected to an inelastic rope. The angular velocity can be selected to obtain a given horizontal speed of the ornithopter in the platform. Finally, the launcher system is made up of two cams to guarantee horizontal stability and a rail to ensure vertical stability. Once the launcher platform is mounted in the desired position, a couple of servo motors moves the cams opening the space between them, so that the ornithopter is eventually released in the air.

C. Validation
Before launching the vehicle, it is important to note that its mass distribution has a crucial effect on the flight performance through the parameter H [see solutions in (8)- (12)], so the ornithopter was consequently calibrated to maintain the same value in all the experiments, preserving its flight stability. An example of the mass distribution is shown in Fig. 5. The experiments were performed outdoors because the ornithopter needs a minimum height and distance to reach the steady-state conditions. The coordinates were 37 • 29'08.2"N 5 • 38'01.5"W, in a sunny day, with 36 • C maximum temperature and 20 km/h maximum wind velocity. The vehicle was successively launched from the crag of 18 meters high in Fig. 6. After the experiments were performed, the position is triangulated from the three points of view, as shown in Fig. 6. This tracking system allows us not to include a GPS, reducing then the mass of the system and improving the gliding conditions. Then, a bundle adjustment is applied to optimize and refine the platform trajectory. First, the velocity is derived from the tracked trajectory. Second, the ornithopter orientation is processed from the IMU with an attitude and heading reference system (AHRS) [20] using an Extended Kalman Filter (EKF). Third, the signal noise is treated with a Butterworth filter.

A. Experimental validation
A set of five similar experiments was performed with the ornithopter of Fig. 2 in gliding configuration at Re = U c c/ν ∼ O(10 5 ), based on the mean chord length and the stade-state velocity, being suitable this inertial dominance to consider the flow as potential. Indeed, the Reynolds number is similar to that of medium-sized birds such as gulls and vultures. The initial conditions prescribed for the experiments, and of course for the simulations, were U b0 0.5, θ 0 95 • , ω 0 = 0 and γ 0 6.5 • , with the wings and tail positions α 0 = 90 • and α 0t = 70 • , respectively, in all the experiments.
For validation purpose, the altitude Z, the speed U b and the pitch θ and gliding γ angles are compared against the numerical results as shown in Fig. 7. The temporal integration of (2)-(5) was performed in Matlab, with a forthorder Runge-Kutta method and a fixed time step dt = 0.1. This time step was selected for providing stable and accurate results compared with other values.  Fig. 7 that the theoretical model captures quite well the trend of the flight magnitudes when compared with the mean of the complete set of experiments. For instance, the temporal evolution of both the altitude Z and the velocity magnitude U b , practically matches the numerical solutions.

Note in
On the other hand, the temporal evolution of the mean experimental values of θ is found to provide less accurate results than other measured variables. This discrepancy is probably due to the accumulative error of the internal state estimator which has been built on the top of a low-cost IMU and runs in a low-power embedded computer. Thus, the integration over time of the inertial measurements leads to the error observed in the figure. Additionally, the pitch angle measurement is noisy and it has been filtered in an AHRS algorithm with an EKF and a low-pass filter, appearing here another source of error.
In any case, the experimental curve roughly follows the reference. Finally, regarding the gliding angle evolution γ, one can notice that the model captures quite well the values obtained experimentally, especially during the last moments, when the steady state is being reached. The width of the shaded zones, indicating the deviation of the experimental mean values, is thought to be due to the gust disturbances appearing on the vehicle during the flight, which sometimes made the trajectory to not to be fully rectilinear.

B. Realistic simulations
Large-scale environments are a key factor in simulating ornithopter models since these have large autonomy and can extend great distances. The maneuverability of these aerial systems require large open spaces to perform turns in a safe way. Commonly used simulators in multi-rotor devices like Gazebo [21] can not recreate physically and visually complex scenes without suffering on performance.
In this work, a sandbox framework was developed with different environments, to allow developers to easily use it and interact with the ornithopter model in Sec. II to manipulate and visualize the flying behavior of the UAV.
This environment has been built on Unreal Engine 4 (UE4), providing powerful rendering and a physics engine that makes it perfectly suitable for robotics simulations. A novel vehicle model has been developed in the Airsim simulator in the UE4 system. The framework architecture in Fig. 8 depicts the implementation of the ornithopter model in the Airsim [22] scheme. The physics engine is divided in two states. If the robot collides with another object, the Airsim physics engine computes the collision response as the Coulomb friction. In any other case, the ornithopter model controls the dynamic behavior. The Airsim sensor model collects physics magnitudes from the environment model providing gravity, magnetic field, air pressure, and density. The controller module is fed with a barometer, gyroscope, accelerometer, and magnetometer. Then a control signal is sent to the ornithopter to close the loop. The visual model developed and the environment are rendered at a constant 60 Hz refresh rate. Therefore, the ornithopter model in (2)-(5) was implemented within the Airsim framework and integrated numerically in time with a fourth-order Runge-Kutta method. Flapping-flight episodes are also explored in this novel tool, computing the thrust force in (13) with the Hankel functions provided by the Boost library.
Simulations were tested in a laptop with an Intel Core i7-7700HQ CPU and NVIDIA GeForce GTX 1070 GPU. Figure 9 shows some snapshots of the simulator with the flying platform in realistic environments, in both gliding and flapping-wing flight configuration.

V. CONCLUSION
This work provides a simple model derived from the linearized potential theory that can be used to understand the aerodynamic behavior of a flapping-wing UAV in gliding and low-amplitude flapping-wing flight configurations. It is shown that the results obtained experimentally agree quite well with those of the model, capturing the overall trend in the flight variables U b , θ and γ, and the altitude Z. Despite the environmental disturbances, the agreement between numerical and experimental results show that the theoretical model could be useful to estimate the flight features of an ornithopter. This fact has direct implications in UAVs because it is useful to know the vehicle dynamics prior to incorporate appropriate control systems to make the flight stable.
On the other hand, the proposed model suffers from certain limitations. First, it is worth to note that the aerodynamic forces involved in the gliding model are only valid for small angles of attack (for wings and tail), corresponding it to practically horizontal flights. Second, if wing strokes are enabled, only small amplitudes are likely to be modeled, guaranteeing no flow separation. Finally, the present model is fully longitudinal, so only the dynamics contained in the vertical plane are addressed.
Concerning the launcher platform, it has taken a significant role during the experimental validation. It was crucial for setting the gliding conditions during experiments. Nevertheless, weather conditions delimited the flight day. The data is publicly available and accessible on the following link https://doi.org/10.5281/zenodo.3408565 for the community.
In addition to the experimental validation, a novel platform incorporating the model equations was implemented within the Airsim framework. This contribution results to be highly efficientproviding the following: • High-frequency calculations to operate real-time hardware in the loop (HITL) simulations and compare the model with real-world flights. • UE4 brings photo-realistic visual rendering which allows collecting a large amount of annotated data with Airsim tools. Consequently, the flexible environment of UE4 with his large online marketplace of realistic models grants a perfect system to gather valuable datasets in a variety of conditions. • The simulator supports both flapping and gliding transitions. As future works, one of the most challenging and immediate tasks would be the experimental validation of the present model switching between gliding flight configuration to flapping-wing episodes. It has quite important implications in the study of the efficiency of the forthcoming flappingwing aerial vehicles. The incorporation of an appropriate control system for this flight configuration will be also needed to prescribe a desirable path to be followed by the ornithoper. It is also of importance the extension of the model to a three-dimensional space in order to perform turning manoeuvres.