Simplified Model for Forward-Flight Transitions of a Bio-Inspired Unmanned Aerial Vehicle

Citation: Sanchez-Laulhe, E.; Fernandez-Feria, R.; Ollero, A. Simplified Model for Forward-Flight Transitions of a Bio-Inspired Unmanned Aerial Vehicle. Aerospace 2022,9, 617. https://doi.org/ 10.3390/aerospace9100617 Academic Editor: Giuseppe Pezzella Received: 14 September 2022 Accepted: 14 October 2022 Published: 18 October 2022 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). aerospace Article Simplified Model for Forward-Flight Transitions of a Bio-Inspired Unmanned Aerial Vehicle Ernesto Sanchez-Laulhe 1,*, Ramon Fernandez-Feria 1and Anibal Ollero 2 1Fluid Mechanics Group, University of Malaga, Dr Ortiz Ramos S/N, 29071 Malaga, Spain 2GRVC Robotics Laboratory, University of Seville, Camino de los Descubrimientos S/N, 41092 Seville, Spain *Correspondence: [email protected] Abstract: A new forward-flight model for bird-like ornithopters is presented. The flight dynamics model uses results from potential, unsteady aerodynamics to characterize the forces generated by the flapping wings, including the effects of the dynamic variables on the aerodynamic formulation. Numerical results of the model, which are validated with flapping flight experimental data of an ornithopter prototype, show that state variables such as the pitch angle and the angle of attack oscillate with the flapping frequency, while their mean values converge towards steady-state values. The theoretical analysis of the system shows a clear separation of timescales between flapping oscillations and transient convergence towards the final forward-flight state, which is used to substantially simplify both the interpretation and the solution of the dynamic equations. Particularly, the asymptotic separation into three timescales allows for dividing the problem into a much simpler set of linear equations. The theoretical approximation, which fits the numerical results, provides a direct look into the influence of the design and control parameters using fewer computational resources. Therefore, this model provides a useful tool for the design, navigation and trajectory planning and control of flapping wing UAVs. Keywords: unsteady aerodynamics; flight dynamics model; bio-inspired UAV; perturbation method; ornithopter 1. Introduction Research on flapping wing Unmanned Aerial Vehicles (UAVs) is becoming of increasing interest due to their potential to improve the efficiency and safety of multi-rotors without the loss of maneuverability of fixed-wing aircrafts [ 1 ]. The ability to glide for long distances gives bird-like ornithopters great energy efficiency, while the flapping wings provide the needed thrust to gain altitude. Such advantages have led to the development of several prototypes of flapping wing bio-inspired UAVs. For instance, RoboRaven [ 2 ] was the first demonstration of a bird-inspired platform performing outdoor aerobatics using independently actuated and controlled wings. Several ornithopter prototypes were developed in [ 3 ], evaluating lift and thrust generation characteristics of different wing designs. The RoBird project [ 4 ] achieved the speed range of predatory birds with similar weights. Currently, the GRIFFIN project is focused on the research and development of bio-inspired UAVs, with the E-flap prototype already presenting an ability to carry a high payload [5]. Predicting and controlling the behavior of a flapping wing UAV constitutes a great challenge. Fully empirical models may be useful [6], but are restricted to the specific characteristics of a particular ornithopter. General models have to be based on the flight physics of the ornithopter, integrating aerodynamic forces and moments generated by the wings into the dynamics of the entire ornithopter. There are already different formulations, most of them focused on very small or Micro Aerial Flapping Wing Vehicles (MAFWVs) [7–9] . Such models combine estimations of the aerodynamic forces for high angles of attack, Aerospace 2022,9, 617. https://doi.org/10.3390/aerospace9100617 https://www.mdpi.com/journal/aerospace Aerospace 2022,9, 617 2 of 22 acting in aerodynamic stalls, with complex flapping kinematics similar to those produced by an insect. From their wing kinematics, MAFWVs obtain the incoming velocity needed for the generation of lift and thrust. For larger UAVs, mechanical limitations do not allow this flight mode, so they need to start from a forward flight motion, e.g., from gliding, to produce the aerodynamic forces. Unsteady aerodynamics based on inviscid (or potential) theory is useful to model these forward flapping-flight regimes. Flapping wing aerodynamics in forward flight has been a research topic since Theodorsen [ 10 ] and Garrick [ 11 ] formulated the lift and thrust for a two-dimensional flapping airfoil within the linear potential theory framework. Theodorsen’s formulation was used by DeLaurier [ 12 ] to model flapping-wing flight but also considered large amplitude wing motion to develop the Modified Strip Theory (MST). Kim et al. [ 13 ] further improved this model to include relatively large angles of attack and dynamic stall effects. Recently, Garrick’s linear inviscid formulation for the thrust force was improved using the vortical impulse theory [ 14 ], which considers the complete unsteady vorticity distribution on the airfoil. Some recent studies have adapted these theoretical two-dimensional formulations to real finite wings [15,16]. Different flight models may integrate these aerodynamic forces into the dynamics of the entire ornithopter. Rigid body approximation was used by Dietl and Garcia [ 17 , 18 ], using a simplification of the modified strip theory. Paranjape et al. [ 19 ] adapted insect-like models to forward flights with cycle-averaged values of the forces. Similar empirical models to compute the aerodynamic forces were integrated into the longitudinal flight dynamics by Taylor and Zhikowski [20] and by Gim et al. [21]. With high amplitude flapping, the rigid body assumption loses validity. Some works have studied multi-body models, considering the interaction between moving and rigid parts [ 22 ]. Most of the multi-body models in the literature do not use theoretical formulations for the aerodynamic forces but instead use empirical approximations [ 6 , 23 ]. However, there are also some works that use modified strip theory for flexible wings [ 24 ], some considering the fluid–structure interaction [ 25 ] but simplified for computational purposes. One of the most relevant current challenges is to develop general flight models for actual ornithopters but simple enough to be used to control the UAV with on-board computers in real-time. To be general enough, such models must be developed from physical theories. Within this general goal, previous work by the authors on the GRIFFIN project has focused on characterizing gliding states and their stability [ 26 , 27 ]. The aim of the present work is to characterize the ornithopter flight transitions from steady gliding states to long-time, or permanent, flapping-flight states after a flapping regime of small amplitude is started. The model may also be applied to the transitions between two different flapping regimes. The work is limited to forward flight in the longitudinal plane, as only symmetrical behavior is studied. Lateral stability issues are not addressed here. This paper presents a rather simple model integrating low-amplitude flapping aerodynamics in the longitudinal dynamics of a bird-inspired UAV. Formulations of the lift [ 10 ] and thrust [ 14 ] from linear potential theory, corrected to account for finite-wing effects, provide a precise estimation of the aerodynamic force generation by the flapping wings, which has been validated against experimental data [ 28 – 31 ]. These forces are included in the Newton–Euler dynamic equations, coupling the vehicle dynamics and the wing aerodynamics. Since the aerodynamic model used is for low-amplitude flapping, the use of rigid body equations is justified. Furthermore, the reduced weight of the wings, normally required by mechanical reasons, minimize the impact of their relative movement in the dynamics of the vehicle. Aerodynamic formulation allows consideration of the unsteady damping for both the wing and the tail, an aerodynamic effect that is relevant in the evolution of the system towards its final flapping-flight state. The dynamics equations are written in the velocity frame to simplify the formulation. These assumptions provide a model that we believe is quite general for a wide range of ornithopters in longitudinal forward flight. Aerospace 2022,9, 617 3 of 22 Despite the simplifications, the numerical integration of the resulting dynamic equations is still computationally expensive. This is due to the quite different timescales; for flapping, the timescale is much smaller than that for the transient phase between, for instance, an initial gliding state and a final flapping-flight state. Therefore, a huge number of time steps is required to obtain accurate numerical solutions in both timescales. Hence, onboard real-time numerical implementation of these equations is not possible. The option of using time-averaged expressions for the lift and thrust is discarded because of the nonlinear coupling of the dynamics and aerodynamics of the system, as flapping oscillations affect the temporal evolution of the time-average quantities and, consequently, the dynamics of the ornithopter towards its final forward flight state. The path followed here is to take advantage of the disparate timescales to develop a multiple scales perturbation method to solve the dynamic equations, as applied, for instance, in [ 32 ] for a much simpler fish swimming problem. For that purpose, we use the dimensionless flapping amplitude (scaledwith the wing’s chord length) asthe small parameter. In summary, the main contribution of this paper is the development of a forward flight model with linear potential aerodynamic results [ 10 , 14 ] for flapping wing force generation and its rigorous simplification based on the separation of the different timescales, providing a useful analytical tool for the design, real-time control and guidance of bird-inspired UAVs. The paper is structured as follows. Section 2presents the physics of the ornithopter flight, with both the Newton–Euler equations used for the dynamics and a summary of the unsteady aerodynamic forces generated by the flapping wings. The details of the aerodynamic forces are given in Appendix A. Section 3shows some numerical results and the experimental validation of the model. Section 4explains the approach followed to obtain the simplified solution of the model and presents the comparison between this simplified analytical solution and the numerical simulations. Finally, in Section 5, the main conclusions are summarized. 2. Flight Model This section presents the dynamic equations and the aerodynamic model used for the ornithopter in flapping flight. The implementation of Newton–Euler equations for a rigid flying vehicle in longitudinal flight is based on previous works [ 26 , 27 , 33 ] and summarized in non-dimensional form in Section 2.1. Although these past works were focused on gliding flight, the models can be generalized for the flapping case by adding the appropriate aerodynamic forces. We use an aerodynamic model for flapping wings in forward flight based on linearized potential aerodynamics. The formulation is given in Section 2.2, together with the aerodynamic model for the fixed but adjustable tail, whose deflection angle δt (see Figure 1) is used as the main flight control parameter. In particular, we shall assume a heaving motion of the wings of the form ˜ h(t) = <h˜ h0eiω˜ ti, (1) where ω and ˜ h0 are the flapping frequency and amplitude, respectively, ˜ t is the time, and < means a real part. A tilde ˜ is used on dimensional quantities for which the same symbol will be used later for their dimensionless counterparts. Although the amplitude of the wing flapping motion varies from the root to the tip, to adjust the formulation, a characteristic amplitude of the flapping wing is defined as that of the chord placed at 1 / 3 of the wing semi-span from the wing tip, as it has been proved to be adequate in previous works [34]. Aerospace 2022,9, 617 4 of 22 Figure 1. Scheme of forces, axis, and control and state variables. All vectors are on the plane ( X , Z ) of the longitudinal flight. 2.1. Non-Dimensional Newton–Euler Equations A sketch of the ornithopter with the different reference frames and the main (dimensional) parameters and state variables are provided in Figure 1. To simplify the problem, we use non-dimensional magnitudes scaled with the following characteristic velocity, length and time: Uc=rmg πρSδt,Lc=c 2,tc=1 ω, (2) where Uc and tc are redefined from [ 27 ] to obtain a non-dimensional velocity of order unity in usual flapping conditions for a given tail deflection δt , and a non-dimensional time in the order unity for flapping oscillations. Note that this formulation is only valid for the flapping mode of the ornithopter, as the flapping frequency is involved in the scaling of the dimensionless time. During a longitudinal flight, one needs three dynamic equations: two for the linear momentum on the plane of flight ( X - Z , see Figure 1) and one for the angular momentum around the Y -axis perpendicular to this plane. They are written as [ 27 ], with some modifications Mk0˙ U=U2(C∗ T−C∗ D−Li∗−ΛC∗ Dt)−δtsin (γ),(3) Mk0U˙ γ=U2(C∗ L+ΛC∗ Lt)−δtcos (γ),(4) (Mk0)2¨ θ=χM2U2{lw[C∗ Lcos (α)−(C∗ T−C∗ D)sin (α)]+ltΛ[C∗ Lt cos (α) + C∗ Dt sin (α)] −hw[C∗ Lsin (α) + (C∗ T−C∗ D)cos (α)]},(5) and the wing’s angle of attack is computed as α=θ−γ . Equations (3) – (5) are written in the trajectory frame, with xT in the direction of the velocity. The non-dimensional parameters appearing in Equations (3)–(5) that differ from those used in [27] are M=2m πρSc,χ=πρSc 23 2Iy,k0=ωc 2Uc, (6) The (modified) force coefficients C∗ L,C∗ T,C∗ D, modeled in Section 2.2 are defined as C∗ L=L πρU2 bc,C∗ T=T πρU2 bc,C∗ D=D πρU2 bc, (7) applied at the wing aerodynamic center. The Lighthill number Li∗=Sb SC∗ Db , models the body drag, applied at the center of gravity. Small flapping amplitude is considered so that the wing aerodynamic center is assumed to be at a constant position ( lw , hw ) in the Aerospace 2022,9, 617 5 of 22 body frame. The force coefficients are dubbed modified because they are not defined in the traditional form but divided by 2 π to better adjust the numerical order of magnitude of the different terms in the equations, as described in Section 4. 2.2. Aerodynamic Characterization For the aerodynamic formulation, the movement of the wing is characterized by the non-dimensional flapping amplitude h0 and the reduced (dimensionless) frequency, k=k0 U , which is affected by the dynamics through the non-dimensional velocity U in relation to the constantly reduced frequency k0, defined in (6). Since we assume small flapping amplitudes, we use the well-known results from twodimensional, linear potential theory, recently modified for the thrust coefficient and also corrected to account for three-dimensional effects. Thus, we use Theodorsen’s lift [ 10 ] for an oscillating, heaving-only two-dimensional airfoil but assume an angle of attack α6= 0 that may vary with time. In fact, all terms coming from changes in dynamic state variables are also considered as they have important effects on the stability of fixed-wing flight [ 27 ] and are also expected to affect the flapping flight. To also account for three-dimensional effects, we follow recent works that separate the finite wing effects on the added mass terms and the circulatory terms of the lift force [ 15 ]. Circulatory terms are affected by the same factor used in fixed wings according to Prandtl’s lifting-line theory, namely A/(A+ 2 ) , where A=b2/S is the aspect ratio of the wing. On the other hand, Smits [ 16 ] proposes a discontinuous factor for the added mass terms, with A/ 2 for A< 2 and just unity for A≥ 2. This last factor is the appropriate one for the wings considered here. We can write the (modified) lift coefficient as a function of the state variables: C∗ L=CLαα+CLheit +¯ CLhe−it Uh0+CL˙ α U˙ α+CL˙ U U2α˙ U+CL˙ θ U˙ θ, (8) where we see the harmonic function related to the heaving motion multiplying h0 . The different coefficients CL... in Equation (8) (the dots refer to the different subscripts) are obtained from Theodorsen’s theory, and they are defined in Appendix A. The formulation of the thrust for an oscillating airfoil in the linear potential limit was first obtained by Garrick [ 11 ]. The theoretical development has been recently corrected in [ 14 ] using the vortex impulse theory. This formulation for the thrust is used here. Finite wing effects from [15,16] are considered, resulting in C∗ T=CTc U2h2 0+CTαα2+CThαeit +¯ CThαe−it Uh0α+CThe2it +¯ CThe−2it U2h2 0+CTh˙ αeit +¯ CTh˙ αe−it U2h0˙ α +CT˙ θ˙ α U2˙ θ˙ α+CT˙ αα U˙ αα +CTh˙ θeit +¯ CTh˙ θe−it U2h0˙ θ+CT˙ θα U˙ θα +CT˙ θ U2˙ θ2+CT˙ Uα U2α2˙ U. (9) Note that the variation of the dynamic variables produces several cross terms. For an isolated wing in a wind tunnel with a heaving movement, just the first four terms of Equation (9) would be present. All the coefficients CT... (dots refer to the different subscripts in Equation (9) ) are defined by the development of the formulation in [ 14 ], and they are also defined in Appendix A. The second lifting surface of the ornithopter, the tail, provides stability and pitch control. To formulate the forces generated by this second surface, delta-wing theory is used [ 35 ]. Bio-inspired UAVs usually have a triangular tail with a reduced aspect ratio similar to those of birds, for which this theory has proven to be appropriate [ 36 , 37 ]. Consequently, the tail’s lift is modeled as [27] C∗ Lt =CLtα(α−δt)+CLt ˙ θ U˙ θ+CLt˙ α U˙ α=At 4(α−δt)−At 4 k0lt U˙ θ+3At 8 k0 U˙ α, (10) Aerospace 2022,9, 617 6 of 22 Finally, drag from all parts of the vehicle has to be considered. Body drag, represented by the modified Lighthill number Li∗ , is considered constant for the normal range of angles of attack. Meanwhile, for the wing and tail, drag is computed as the contributions of a friction drag, which is also considered constant in the range of applicability of the model, and an induced drag is associated with the corresponding lift coefficients C∗ L and C∗ Lt , which is computed from finite wing theory. In summary, C∗ D=C∗ D0+C∗ Di,C∗ Di =2C∗2 L A, (11) C∗ Dt =C∗ D0t+C∗ Dit,C∗ Dit =2C∗2 Lt At. (12) 3. Numerical Results and Experimental Validation This section presents, in Section 3.1, some numerical results of the model equations to describe the structure of the solution, which will guide the timescales analysis of the perturbation solution developed in Section 4. Before that, in Section 3.2, these numerical results are validated against experimental data from the flapping flight of an ornithopter prototype developed within the GFIFFIN project [5]. 3.1. Numerical Results To obtain the dynamic evolution from the gliding condition to the final flapping state, the equations described in Section 2are numerically integrated. We use the values of the parameters from the prototype described in [ 5 ], which are shown in Table 1for three flapping frequencies f relevant for bird-scale flight (actually, only k0 depends on f ). Note that they are of order unity, except for the friction drag (Li∗). Table 1. Non-dimensional flight parameters for three flapping frequencies. fMk0MM2χΛlwltΛhwLi∗A 2Hz 0.79 2.54 2.12 0.25 0.55 −1.16 0.38 0.0048 5.14 5Hz 1.98 2.54 2.12 0.25 0.55 −1.16 0.38 0.0048 5.14 7Hz 2.77 2.54 2.12 0.25 0.55 −1.16 0.38 0.0048 5.14 With these parameters, we proceed to numerically integrate the system of differential Equations (3) – (5) . Figure 2shows the evolution of the non-dimensional flight velocity U with both non-dimensional time (a) and dimensional time (b) for the three frequencies in Table 1. The flapping amplitude considered is h0= 0.1. To better compare the three cases, velocities have been normalized with their final values Uf , which are 0.9136, 1.0444 and 1.0792, respectively, and are all of order unity, as mentioned in Section 2.1. From Figure 2a, it is observed that the timescale of the transient phase increases with the frequency, while Figure 2b shows that the transient velocity is just slightly affected by the frequency when the dimensional time is used. Thus, though the timescale related to flapping oscillations is constant as a consequence of the non-dimensionalization of time with the flapping frequency, the timescale related to the transition between states varies with the frequency for the same reason. Aerospace 2022,9, 617 7 of 22 0 100 200 300 400 500 600 700 800 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 2 Hz 5 Hz 7 Hz (a) 0 10 20 30 40 50 60 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 2 Hz 5 Hz 7 Hz (b) Figure 2. Transient phase of non-dimensional velocity U for three different frequencies. ( a ) Normalized velocity vs. non-dimensional time; (b) normalized velocity vs. time in seconds. Figure 3shows the evolution of the angular variables θ and α with time for a frequency of 5Hz . Notice that the pitch angle amplitudes in the transient oscillations are significantly higher than those of the angle of attack. However, the oscillations due to the flapping motion are both quite similar. All the amplitudes, both in the transient oscillation and during each flapping cycle, are small enough to be in the adequate range for the aerodynamic models during the simulation [ 28 – 31 ]. Notice also that the aerodynamic models have no limitation in relation to the frequencies. 0 200 400 600 800 1000 1200 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 (a) Evolution of the pitch angle 0 200 400 600 800 1000 1200 0 0.05 0.1 0.15 0.2 0.25 (b) Evolution of the angle of attack Figure 3. Evolution of θ(a) and α(b) with non−dimensional time. Frequency 5Hz. Numerical integration of the dynamic equations provides the most accurate tool to obtain the evolution of the state variables. However, state variables (velocity, pitch angle and angle of attack) oscillate with the flapping frequency but also evolve much more slowly in time during the transient phase. This disparity of timescales forces the numerical integration to be developed with a very small time step to accurately catch these oscillations, thus requiring a huge number of time steps to solve the slow evolution. Therefore, simulations of the entire transient phase have an excessive computational cost, which makes on-board real-time computations infeasible. However, this same disadvantage associated with the different timescales can be used to develop simplifications of the model. Consequently, this paper proposes a perturbation approach (Section 4), separating scales with the flapping amplitude h0 as the small parameter. The goal is to provide a system of equations that can be computed in real time by on-board computers. 3.2. Experimental Validation Experiments with the ornithopter prototype (described in [ 5 ]; see Figure 4) have been carried out in a 20 m × 15 m tracking zone inside the GRVC Robotics Laboratory, with five Aerospace 2022,9, 617 8 of 22 markers on the body for the cameras to track the flight. The limited space of the flight area is inconvenient for achieving a long flight time, but its location inside a building without wind helps for the repeatability and reliability of the results. The lab is also equipped with a launcher to fix the initial conditions (also shown in Figure 4). Figure 4. Ornithopter prototype used in the experiments. The prototype has a flapping amplitude h0≃ 0.8 (at 1 / 3 of the semi-span from the wing tip), somewhat out of the range of the linear aerodynamics models. However, it has been shown that these models provide a good approximation of the aerodynamic forces even for these flapping amplitudes [ 14 , 28 – 31 ] (see also the present comparison with experimental data reported below in this section). Although the model does not take into account wing flexibility, the wing used here is more rigid than that reported in [ 5 ], and the effect of flexibility on the forces generated by the flapping wing is small for the flapping frequencies considered, which are much smaller than the characteristic frequency of the structure. Tail flexibility has not been considered either. In addition, due to the characteristics of the prototype’s tail, the coefficient CLt ˙ θis very small and is neglected. There are three controls on the ornithopter. The vertical tail has a simple control to follow a straight line in the projection on the horizontal plane. We follow the diagonal of the flight area to maximize flight time. To consider the flight in open loop, the horizontal tail remains in a fixed position by sending a constant signal to the servomotor. The third control of the ornithopter generates the flapping frequency, which also receives a constant signal. However, this signal generates a constant current, which does not guarantee that the flapping motion is exactly harmonic, though the error in relation to a harmonic motion is shown to remain very small. The flapping frequency is not known a priori, but it can be extracted from the oscillations in the measured variables. We consider the flight from the moment the flight starts. Due to electronic delays, the ornithopter starts to fly between 0.5 and 1 m forward from the launcher. The initial conditions change slightly from the launch to that point, so the conditions defined in the launcher are approximate. Figure 5shows that the numerical results of the model for the flight velocity U and the pitch angle θ are very close to the experimental measurements. As we see, the flight time is only around 2s due to the reduced space in the tracking zone. However, results for this short time flight are significant, as at least one transient oscillation is caught, with good agreement between experimental and theoretical frequencies and amplitudes of both fast (flapping) and transient scales. More particularly, the amplitude of the small oscillations with the flapping frequency have the same order for the airspeed, but they are slightly different for the pitch angle. This difference could be explained by the impact of the inertial forces, which would move the aerodynamic center of the wing backward, closer to the center of gravity. On the other hand, the first transient oscillation is very well captured by the model for both Uand θ. Aerospace 2022,9, 617 9 of 22 0 2 4 6 8 10 3.5 4 4.5 5 5.5 6 6.5 7Numerical Simulation Experimental Data (a) 0 2 4 6 8 10 -0.2 0 0.2 0.4 0.6 0.8 1Numerical Simulation Experimental Data (b) Figure 5. Comparison between numerical and experimental evolutions. ( a ) Airspeed; ( b ) pitch angle. 4. Analytical Perturbation Solution and Its Comparison with the Numerical Solution In this section, we present a perturbation solution to the model equations. First, in Section 4.1 , the different timescales and orders of magnitude of the different parameters are analyzed with the help of the above numerical solutions and theoretical considerations. The solution is fully derived in Section 4.2 and then compared to the numerical results of the model equations in Section 4.3. 4.1. Scales Analysis It is clear from the above numerical and experimental results that the problem has different timescales, so it seems appropriate to use a multiple-scale perturbation method, whose different temporal scales will be defined here in terms of the small flapping amplitude, e≡h0 1. Before that, some assumptions have to be made regarding the order of magnitude, in terms of e , of the non-dimensional variables and the remaining parameters for a typical ornithopter, which are enumerated next. • Small flight angles: γ∼θ∼α∼δt∼e . Reduced angles of attack are required by the present aerodynamic formulation, as well as small deflections of the tail. All of them are assumed to be in the same order of magnitude as the flapping amplitude e . Small γ and θ mean that aggressive climbing trajectories are not considered. Moreover, small flapping amplitude implies a limited thrust, which is proportional to e2 (see below), so reaching large flight path angles would be possible only with high reduced frequencies, which are difficult to obtain in bird-scale ornithopters. • Lift coefficients: C∗ L∼C∗ Lt ∼e . This is a consequence of the formulation of the lift coefficients, Equations (8) and (10) , proportional to the aerodynamic angles and the flapping amplitude, as all the aerodynamic parameters are of order unity. • Thrust coefficient: C∗ T∼e2 . Clearly from the Formulation (9) , with products of the aerodynamic angles and the flapping amplitude. • Drag coefficients: C∗ D∼C∗ Dt ∼Li∗∼e2 . This is clear for the induced drag, as it is proportional to the square of the lift coefficient. On the other hand, friction drag cannot be larger than thrust to maintain flight. • Non-dimensional inertia and reduced frequency: M2χ∼ Mk0∼ 1. Given in Table 1. Even though the frequency is a control variable, we can consider both terms of order unity in the typical range of bird-like ornithopters. • All the remaining non-dimensional parameters either do not affect the structure of the solution or are order unity. With these considerations, we can analyze the timescales. Two of them are clearly recognized in Figure 2. The first one, corresponding to the flapping oscillations, is of order unity (due to the definition of the characteristic time tc ). The order of magnitude of the other one, corresponding to the transient phase, can be obtained by analyzing the order of the different terms in Equations (3) and (4) during the transient phase. The asymptotic terms are neglected in this analysis as they are compensated by themselves. Thus, writing Aerospace 2022,9, 617 16 of 22 0 200 400 600 800 1000 1200 1400 0.6 0.8 1 1.2 1.4 Numerical Simulation Analytical Approximation (a) 0 200 400 600 800 1000 1200 1400 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 Numerical Simulation Analytical approximation (b) Figure 7. Comparison between theoretical and numerical transient phase at 5Hz. ( a ) Airspeed evolution; (b) pitch angle evolution. Figure 8shows how the error in the transient approximation is not important for the trajectories, accumulating just a very small error in the large-time longitudinal path. Note how the thrust provided by the wing in this particular case is not enough to climb up, and the ornithopter follows a descending trajectory. To achieve more aggressive trajectories, an additional pitching movement would have to be included with the flapping motion. Most prototypes include this movement, either by an active pitching mechanism or by means of wing flexibility. 0 100 200 300 400 500 600 700 800 900 -80 -70 -60 -50 -40 -30 -20 -10 0 10 Numerical Simulation Analytical Approximation Figure 8. Comparison between the trajectories at 5 Hz. 5. Summary and Conclusions A new physical model for the forward flight of flapping wing UAVs has been presented. The model is restricted to low-amplitude flapping, appropriate for UAV cruising, using results from potential aerodynamic theory. Wholly unsteady aerodynamics models are employed instead of the quasi-steady approximations normally used [ 18 , 25 ], obtaining a more precise characterization of the non-stationary forces generated by the flapping wings. This new model is general for any bird-like ornithopter with a simple flapping (heaving) motion, and its results are validated with experimental data from the longitudinal flight of an ornithopter prototype. The low-amplitude flapping approximation allows simplifications of the dynamic characteristics of the ornithopter. However, the numerical simulations with the model still have a high computational cost due to the small time step needed to accurately capture the flapping oscillations along very long transient periods. The resources needed for realtime on-board computations would exceed the limited payload of ornithopters. Hence, this paper proposes a simplification using multiple-scale perturbation methods in terms of the small flapping amplitude. The solution consists of a transient oscillatory phase, which converges to an asymptotic oscillatory regime with a much smaller amplitude and Aerospace 2022,9, 617 17 of 22 a much larger flapping frequency. Thus, the perturbation solution is built up with three timescales, the shortest one associated with the flapping frequency and an intermediate one corresponding to the transient oscillations towards the long-time behavior, which is reached in a much slower timescale. The perturbation solution is obtained analytically up to the second and third order in the small parameter, depending on the state variable. In addition to its utility for on-board computations in real-time for trajectory planning and control, the analytical solution presented here is of interest, providing a simple and straightforward characterization of the final flapping-flight cruising state. The validity of the approach is assessed by comparing it with numerical results previously validated against experimental data. It is shown that the large-time asymptotic values are very well caught by the approximate solution for both the main values and the amplitude and phase of the oscillations. The approximate solution for the transient phase proves the validity of the selected timescales, although its agreement with the exact numerical solution is less accurate. However, the final accumulative error in the trajectory is very small, as the larger transient errors disappear as the asymptotic values are reached. The model is limited to bird-scale flapping-wing robots for the assumptions made on the order of magnitude of the variables to be valid. These assumptions are essential for the adequate development of the perturbation method. A different scaling would lead to changes in the perturbation equations. Thus the model may fail for very small robots such as hummingbirds and insect-like UAVs. Another limitation is associated with the aerodynamic theory used, valid for small angles of attack and low-amplitude flapping. However, the models for the aerodynamic forces used here have been previously validated with experimental results, showing good agreement even for not-so-small flapping amplitudes in the present case of heaving motion [14,28–31], and the results of the dynamic equations with these aerodynamic models have been validated in the present work against experimental measurements of the longitudinal flapping flight of an ornithopter prototype. Additionally, a rigid wing with just a flapping (heaving) movement is considered. Pitching, both active and/or passive by flexibility, is usual for bird-like vehicles because, as also shown here, flapping rigid wings produce limited thrust. With such low thrust, ornithopters need a very low drag to maintain sustained flight. Therefore, a combination of pitching and flapping movements is the most interesting path to extend this work. Considering both active pitching mechanisms and passive flexibility pitching, the model would be effective for all actual designs of flapping wing bird-like UAVs. Finally, it is worth mentioning that another interesting path to extend this work would be the inclusion of lateral dynamics. Ornithopters do not usually turn like fixed-wing aircrafts, as the inclusion of anti-symmetrical control on the wing causes high mechanical complexity, especially in bird-sized UAVs. Then, most of them rely on tail configurations to control the lateral trajectory, obtaining different dynamics from those of conventional aircraft. Author Contributions: Conceptualization, E.S.-L., R.F.-F. and A.O.; methodology, E.S.-L. and R.F.-F.; software, E.S.-L.; validation, E.S.-L.; formal analysis, E.S.-L.; investigation, E.S.-L. and R.F.-F.; data curation, E.S.-L.; writing—original draft preparation, E.S.-L.; writing—review and editing, R.F.-F. and A.O.; visualization, E.S.-L.; supervision, R.F.-F. and A.O.; funding acquisition, A.O. All authors have read and agreed to the published version of the manuscript. Funding: The authors acknowledge support from the Advanced Grant of the European Research Council GRIFFIN, Action 788247. E.S.-L. also acknowledges his predoctoral contract at the Universidad de Málaga. Data Availability Statement: The data presented in this study are available on request from the corresponding author. Acknowledgments: Authors acknowledge the support of Mario Hernández in the experimental validation. Conflicts of Interest: The authors declare no conflict of interest. Aerospace 2022,9, 617 18 of 22 Appendix A. Aerodynamic Coefficients in Equations (8)and (9) CLα=A A+2F(k), (A1) CLh=k0 2k 2+A A+2(G(k)−F(k)i), (A2) CL˙ α=CL˙ U=k0 2, (A3) CL˙ θ=−A A+2F(k)k0(lw−1), (A4) where F(k)and G(k)are the real and imaginary parts of Theodorsen’s function [10]. CTc=−k2 0 π G1(k)A A+2, (A5) CTα=−2 π G1(k)A A+2−CLα, (A6) CThα=k0 π(F1(k)+2iG1(k)) A A+2−CLh, (A7) CTh=−ik2 0 2π(F1(k)+iG1(k)) A A+2, (A8) CTh˙ α=ik2 0 4−k2 0A 2A+4(iF2(k)+G2(k)), (A9) CT˙ θ˙ α=k2 0 2lw−1 2−k2 0(lw−1)F2(k)A A+2, (A10) CT˙ αα =−k0 2+k0F2(k)A A+2−CL˙ α, (A11) CTh˙ θ=−k2 0 π iG1(k)A A+2(lw−1)+k2 0 π (−iG1(k)−F1(k))A A+2lw−1 2, (A12) CT˙ θα =2k0 π G1(k)A A+2(lw−1)+2k0 π G1(k)A A+2lw−1 2−CL˙ θ, (A13) CT˙ θ=−2k2 0 π G1(k)A A+2(lw−1)lw−1 2(A14) CT˙ Uα=−CL˙ U, (A15) where F1(k) and G1(k) are the real and imaginary parts of the modified function defined in [14], while F2(k)and G2(k)are given by F2(k) = G(k) + 2 π(kG1(k)−F1(k)),G2(k) = F(k) + 2 π(kF1(k) + G1(k)). (A16) Appendix B. Higher Order Perturbation Terms System of algebraic equations for T3and A3: M1T3 A3=b1(A17) Aerospace 2022,9, 617 19 of 22 being M1(1,1) = iMk0−CL˙ θ−ΛCLt ˙ θ(A18) M1(1,2) = (−i(Mk0+CL˙ α+ΛCLt˙ α)−V0(CLα+ΛCLtα)) (A19) M1(2,1) = − k2 0 χ+iV0lwCL˙ θ+ltΛCLt ˙ θ!(A20) M1(2,2) = −V2 0(lwCLα+ltΛCLtα)+iV0(lwCL˙ α+ltΛCLt˙ α), (A21) and b1(1) = 2V1(CLα+ΛCLtα)A1+V1 V0CLh−iT1Mk0−CL˙ θ−ΛCLt ˙ θ+A1(Mk0+CL˙ α+ΛCLt˙ α)(A22) b1(2) = 2V0V1A1(lwCLα+ltΛCLtα)+V1iT1lwCL˙ θ+ltΛCLt ˙ θ+iA1(lwCL˙ α+ltΛCLt˙ α) +lwV1CLh−V0hwA0V0CTα+CLα2A1+CThα+CLh+(CL˙ α+CT˙ αα )iA1+CL˙ θ+CT˙ θα iT1(A23) Note that T3and A3have a transient term. Its solution can be written as T3=T3s+V1t(τ1,τ2)zT3(A24) A3=A3s+V1t(τ1,τ2)zA3(A25) where zT3and zA3are complex and constant. System of algebraic equations for T4and A4: M2T4 A4=b2(A26) being M2(1,1) = 2iMk0−CL˙ θ−ΛCLt ˙ θ(A27) M2(1,2) = (−2i(Mk0+CL˙ α+ΛCLt˙ α)−V0(CLα+ΛCLtα)) (A28) M2(2,1) = − 4k2 0 χ+2iV0lwCL˙ θ+ltΛCLt ˙ θ!(A29) M2(2,2) = −V2 0(lwCLα+ltΛCLtα)+2iV0(lwCL˙ α+ltΛCLt˙ α), (A30) and b2(1) = 0 (A31) b2(2) = −hwV2 0CTα+CLα+V0A1CThα+CLh+(CL˙ α+CT˙ αα )iA1+CL˙ θ+CT˙ θα iT1+CTh +iCTh˙ αA1−CT˙ θ˙ αA1T1+iCTh˙ θT1−CT˙ θ(A32) Expressions for T2s,A5sand V2s: Aerospace 2022,9, 617 20 of 22 T2s=A2s+2V0V1 δ∗ t(CTα−ΛCDtα)A2 0+2A1¯ A1−ΛCDtαδtδ∗ tA0−C∗ D0−Li∗−ΛC∗ D0t +V1s δ∗ tCThα¯ A1+¯ CThαA1+iCT˙ θα (T1¯ A1−¯ T1A1) +V2 0 δ∗ t(CTα−ΛCDtα)(A0A2s+A1¯ A3s+¯ A1A3s)−ΛCDtαδtδ∗ tA2s +V0 δ∗ tCThα¯ A3s+¯ CThαA3s+iCT˙ θα (T3s¯ A1−¯ T3sA1+T1¯ A3s−¯ T1A3s)+i δ∗ t¯ CTh˙ αA3s−CTh˙ α¯ A3s +CT˙ θ˙ α δ∗ t (T3s¯ A1−¯ T3sA1+T1¯ A3s−¯ T1A3s)+i δ∗ t¯ CTh˙ θT3s−CTh˙ θ¯ T3s+2CTθ δ∗ t (T1¯ T3s+¯ T1T3s)(A33) A5s=−1 V2 0(lwCLα+ltΛCLtα)(2V0V3(lwCLα+ltΛCLtα)¯ A1+2V0¯ V3(lwCLα+ltΛCLtα)A1 +V3lw¯ CLh−i¯ A1(lwCL˙ α+ltΛCLt˙ α)−i¯ T1lwCL˙ θ+ltΛCLt ˙ θ +¯ V3lwCLh−iA1(lwCL˙ α+ltΛCLt˙ α)−iT1lwCL˙ θ+ltΛCLt ˙ θ +2V0V1s(lwCLα+ltΛCLtα)A2s−hw(CTα+CLα)A2 0+2A1¯ A1−ΛCDtαδtδ∗ tA0−C∗ D0 −V1shwCLh+CThα¯ A1+¯ CLh+¯ CThαA1+iCL˙ θ+CT˙ θα (T1¯ A1−¯ T1A1) −hwi¯ CTh˙ αA3s−CTh˙ α¯ A3s+CT˙ θ˙ α(A3s¯ T1+¯ A3sT1+A1¯ T3s+¯ A1T3s)+i¯ CTh˙ θT3s−CTh˙ θ¯ T3s −hwCT˙ θ(T1¯ T3s+¯ T1T3s)+ilwCL˙ U(¯ A1V3+A1¯ V3)+V2 0lw−CLα 2+CTαA3 0+6A0A1¯ A1+A0C∗ D0 +V2 0ltΛ−CLtα 2−CDtαA3 0+6A0A1¯ A1+CLtα 2−CDtαδtA2 0+A1¯ A1δ∗ t+A0C∗ D0 −2hw(CTα+CLα)(A0A2s+A1¯ A3s+¯ A1A3s)) (A34) V2s=−V1s V0+1 2V0(CLαA0+ΛCLtα(A0−δ∗ t))(iMk0(¯ V3(T1−A1)−V3(¯ T1−¯ A1)) −V3¯ CLh−i¯ A1(CL˙ α+ΛCLt˙ α)−i¯ T1CL˙ θ+ΛCLt ˙ θ−2V0V3(CLα+ΛCLtα)¯ A1 −¯ V3CLh−iA1(CL˙ α+ΛCLt˙ α)−iT1CL˙ θ+ΛCLt ˙ θ−2V0¯ V3(CLα+ΛCLtα)A1 −iCL˙ U(¯ A1V3−A1¯ V3)−2V0V1sA2s(CLα+ΛCLtα)−V2 0A5s(CLα+ΛCLtα) −δ∗ tA0 2+A1¯ A1+T0s 2+T1¯ T1−A0T0s−A1¯ T1−¯ A1T1 (A35) Damping of the transient phase ξ ξ=−CL˙ UA0 2ΓMk0(1+∆lw)−ω2 1∆ 2ΓM2χ−CThα¯ A1+¯ CThαA1+iCT˙ θα (T1¯ A1−¯ T1A1)+i¯ CTh˙ αzA3−CTh˙ α¯ zA3 2 −CT˙ θ˙ α(zA3¯ T1+¯ zA3T1+zT3¯ A1+¯ zT3A1)+i¯ CTh˙ θzT3+CTh˙ θ¯ zT3+2CT˙ θ(T1¯ zT3+¯ T1zT3) 2 −V0 2CThα¯ zA3+¯ CThαzA3−iCT˙ θα (−zA3¯ T1+¯ zA3T1−A1¯ zT3+¯ A1zT3)−δ∗ t(T0s−A0) 2Γ −V0(CTα−ΛCDtα)A2 0+2A1¯ A1−ΛCDtαδtδ∗ tA0−C∗ D0−Li∗−ΛC∗ D0t−V0CT˙ θα A0ω2 1 2Mk0 −V2 0 2(CTα−ΛCDtα)(−iω1A0zA2+A1¯ zA3+¯ A1A3s)+iω1ΛCDtαδtδ∗ tzA2+δ∗ tiω1zA2 2(A36) Aerospace 2022,9, 617 21 of 22 and correction of the transient frequency with τ2,ω2 ω2=V1s Γω1(CLαA0+ΛCLtα(A0−δ∗ t)) +CLα(1+∆lw)+ΛCLtα(1+∆lt) ΓV0 ω1A2s−V0V1izA4 +CL˙ θ(1+∆lw)+ΛCLt ˙ θ(1+∆lt)V1sω1 2ΓMk0−∆V0 Γω1(CTα+CLα)A2 0+2A1¯ A1−C∗ D0 −∆hw 2Γω1CLh+CThα¯ A1+¯ CLh+¯ CThαA1+iCL˙ θ+CT˙ θα (T1¯ A1−¯ T1A1)+i¯ CTh˙ αzA3−CTh˙ α¯ zA3 −∆hw 2Γω1CT˙ θ˙ α(zA3¯ T1+¯ zA3T1+zT3¯ A1+¯ zT3A1)+i¯ CTh˙ θzT3+CTh˙ θ¯ zT3+2CT˙ θ(T1¯ zT3+¯ T1zT3) −∆V2 0hw Γω1(CTα+CLα)(−A0iω1zA2+A1¯ zA3+¯ A1zA3)−∆V0hwω1 2ΓMk0CT˙ θα +CL˙ θA1−V1sω1 2Γ −∆V0hw 2Γω1CThα+CLh¯ zA3+¯ CThα+¯ CLhzA3−iCT˙ θα +CL˙ θ(−zA3¯ T1+¯ zA3T1−A1¯ zT3+¯ A1zT3)(A37) with A2(τ1,τ2)defined as A2t(τ1,τ2) = zA2T0t(τ2)eΦT+iω1τ1;zA2=−iω1 Mk0V0 lwCL˙ θ+ΛltCLt ˙ θ lwCLα+ΛltCLtα (A38) and the constant terms Γ=V01−CL˙ θ+ΛCLt ˙ θ Mk0+lwCL˙ θ+ΛltCLt ˙ θ Mk0 ∆(A39) ∆=CLα+ΛCLtα lwCLα+ΛltCLtα (A40) References 1. Mueller, T.J. Fixed and Flapping Wing Aerodynamics for Micro Air Vehicle Applications; AIAA: Notre Dame, IN, USA, 2001. 2. Gerdes, J.; Holness, A.; Perez-Rosado, A.; Roberts, L.; Greisinger, A.; Barnett, E.; Kempny, J.; Lingam, D.; Yeh, C.H.; Bruck, H.A.; et al. Robo Raven: A flapping-wing air vehicle with highly compliant and independently controlled wings. Soft Robot. 2014 , 1, 275–288. [CrossRef] 3. Srigrarom, S.; Chan, W.-L. Ornithopter Type Flapping Wings for Autonomous Micro Air Vehicles. Aerospace 2015 ,2, 235–278. [CrossRef] 4. Folkertsma, G.A.; Straatman, W.; Nijenhuis, N.; Venner, C.H.; Stramigioli, S. Robird: A robotic bird of prey. IEEE Robot. Autom. Mag. 2017,24, 22–29. [CrossRef] 5. Zufferey, R.; Tormo-Barbero, J.; Guzmán, M.M.; Maldonado, F.J.; Sanchez-Laulhe, E.; Grau, P.; Pérez, M.; Acosta, J.Á.; Ollero, A. Design of the High-Payload Flapping Wing Robot E-Flap. IEEE Robot. Autom. Lett. 2021,6, 3097–3104. [CrossRef] 6. Grauer, J.; Ulrich, E.; Hubbard, J., Jr.; Pines, D.; Humbert, J.S. Testing and system identification of an ornithopter in longitudinal flight. J. Aircr. 2011,48, 660–667. [CrossRef] 7. Taha, H.; Hajj, M.R.; Nayfeh, A. Flight dynamics and control of flapping-wing MAVs: A review. Nonlinear Dyn. 2012 ,70, 907–939. [CrossRef] 8. Wang, Z.J. Dissecting insect flight. Annu. Rev. Fluid Mech. 2005,37, 183–210. [CrossRef] 9. Hefler, C.; Kang, C.; Qiu, H.; Shyy, W. Distinct Aerodynamics of Insect-Scale Flight; Cambridge University Press: Cambridge, UK, 2021. [CrossRef] 10. Theodorsen, T. General Theory of Aerodynamic Instability and the Mechanism of Flutter; NACA Technical Report TR 496; NACA: Washington, DC, 1935. 11. Garrick, I.E. Propulsion of a Flapping and Oscillating Airfoil; NACA Technical Report TR 567; NACA: Washington, DC, 1936. 12. DeLaurier, J.D. An aerodynamic model for flapping-wing flight. Aeronaut. J. 1993,97, 125–130. 13. Kim, D.K.; Lee, J.S.; Lee, J.Y.; Han, J.H. An aeroelastic analysis of a flexible flapping wing using modified strip theory. In Active and Passive Smart Structures and Integrated Systems 2008; International Society for Optics and Photonics: SPIE 2008; Volume 6928, 477-484. 14. Fernandez-Feria, R. Linearized propulsion theory of flapping airfoils revisited. Phys. Rev. Fluids 2016,1, 084502. [CrossRef] 15. Ayancik, F.; Zhong, Q.; Quinn, D.B.; Brandes, A.; Bart-Smith, H.; Moored, K.W. Scaling laws for the propulsive performance of three-dimensional pitching propulsors. J. Fluid Mech. 2019,871, 1117–1138. [CrossRef] 16. Smits, A.J. Undulatory and oscillatory swimming. J. Fluid Mech. 2019,874, 1–44. [CrossRef] Aerospace 2022,9, 617 22 of 22 17. Dietl, J.M.; Garcia, E. Stability in ornithopter longitudinal flight dynamics. J. Guid. Control. Dyn. 2008 ,31, 1157–1163. [CrossRef] 18. Dietl, J.M.; Garcia, E. Ornithopter optimal trajectory control. Aerosp. Sci. Technol. 2013,26, 192–199. [CrossRef] 19. Paranjape, A.A.; Dorothy, M.R.; Chung, S.J.; Lee, K.D. A flight mechanics-centric review of bird-scale flapping flight. Int. J. Aeronaut. Space Sci. 2012,13, 267–281. [CrossRef] 20. Taylor, G.K.; ˙ Zbikowski, R. Nonlinear time-periodic models of the longitudinal flight dynamics of desert locusts Schistocerca gregaria. J. R. Soc. Interface 2005,2, 197–221. [CrossRef] [PubMed] 21. Gim, H.; Lee, B.; Huh, J.; Kim, S.; Suk, J. Longitudinal system identification of an avian-type UAV considering characteristics of actuator. Int. J. Aeronaut. Space Sci. 2018,19, 1017–1026. [CrossRef] 22. Grauer, J.A.; Hubbard, J.E., Jr. Multibody model of an ornithopter. J. Guid. Control. Dyn. 2009 ,32, 1675–1679. doi: 10.2514/1.43177. [CrossRef] 23. Amini, M.A.; Ayati, M.; Mahjoob, M. A Simplified Model, Dynamic Analysis and Force Estimation for a Large-scale Orinthopter in Forward Flight Based on Flight Data. J. Bionic Eng. 2020,17, 989–1008. [CrossRef] 24. Pfeiffer, A.T.; Lee, J.S.; Han, J.H.; Baier, H. Ornithopter flight simulation based on flexible multi-body dynamics. J. Bionic Eng. 2010,7, 102–111. [CrossRef] 25. Lee, J.S.; Kim, J.K.; Kim, D.K.; Han, J.H. Longitudinal flight dynamics of bioinspired ornithopter considering fluid-structure interaction. J. Guid. Control. Dyn. 2011,34, 667–677. [CrossRef] 26. Martín-Alcántara, A.; Grau, P.; Fernandez-Feria, R.; Ollero, A. A Simple Model for Gliding and Low-Amplitude Flapping Flight of a Bio-Inspired UAV. In Proceedings of the 2019 International Conference on Unmanned Aircraft Systems (ICUAS), Atlanta, GA, USA, 11–14 June 2019; pp. 729–737. [CrossRef] 27. Sanchez-Laulhe, E.; Fernandez-Feria, R.; Acosta, J.Á.; Ollero, A. Effects of Unsteady Aerodynamics on Gliding Stability of a Bio-Inspired UAV. In Proceedings of the 2020 International Conference on Unmanned Aircraft Systems (ICUAS), Athens, Greece, 1–4 September 2020; pp. 1596–1604. [CrossRef] 28. Ol, M.V.; Bernal, L.; Kang, C.-K.; Shyy, W. Shallow and deep dynamic stall for flapping low Reynolds number airfoils. Exp. Fluids 2009,46, 883–901. [CrossRef] 29. Baik, Y.S.; Bernal, L.P.; Granlund, K.; Ol, M. V. Unsteady force generation and vortex dynamics of pitching and plunging aerofoils. J. Fluid Mech. 2012,709, 37–68. [CrossRef] 30. Fernandez-Feria, R. Note on optimum propulsion of heaving and pitching airfoils from linear potential theory. J. Fluid Mech. 2017,826, 781–796. [CrossRef] 31. Alaminos-Quesada, J. Limit of the two-dimensional linear potential theories on the propulsion of a flapping airfoil in forward flight in terms of the Reynolds and Strouhal number. Phys. Rev. Fluids 2021,6, 123101. [CrossRef] 32. Sánchez-Rodríguez, J.; Raufaste, C.; Argentina, M. A minimal model of self propelled locomotion. J. Fluids Struct. 2020 ,97, 103071. [CrossRef] 33. Lopez-Lopez, R.; Perez-Sanchez, V.; Ramon-Soria, P.; Martin-Alcantara, A.; Fernandez-Feria, R.; Arrue, B.C.; Ollero, A. A Linearized Model for an Ornithopter in Gliding Flight: Experiments and Simulations. In Proceedings of the 2020 IEEE International Conference on Robotics and Automation (ICRA), Paris, France, 2020; pp. 7008–7014. 34. Wang, Z.J.; Birch, J.M.; Dickinson, M.H. Unsteady forces and flows in low Reynolds number hovering flight: two-dimensional computations vs robotic wing experiments. J. Exp. Biol. 2004,207, 449–460. [CrossRef] 35. Jones, R.T. Properties of Low-Aspect-Ratio Pointed Wings at Speeds below and above the Speed of Sound; NACA Technical Report TR 835; NACA: Washington, DC, 1946. 36. Thomas, A.L.R. On the aerodynamics of birds’ tails. Philos. Trans. R. Soc. Lond. Ser. B Biol. Sci. 1993,340, 361–380. [CrossRef] 37. Evans, M.R. Birds’ tails do act like delta wings but delta-wing theory does not always predict the forces they generate. Proc. R. Soc. Lond. Ser. B Biol. Sci. 2003,270, 1379–1385. [CrossRef] 38. Kevorkian, J.; Cole, J.D. Perturbation Methods in Applied Mathematics; Springer: New York, NY, USA, 1981. 39. Ruiz, C.; Acosta, J.; Ollero, A. Aerodynamic reduced-order Volterra model of an ornithopter under high-amplitude flapping. Aerosp. Sci. Technol. 2022,121, 107331. [CrossRef]