Next Article in Journal
A New Method for Analyzing the Aero-Optical Effects of Hypersonic Vehicles Based on a Microscopic Mechanism
Next Article in Special Issue
Lagrange Optimization of Shock Waves for Two-Dimensional Hypersonic Inlet with Geometric Constraints
Previous Article in Journal
Fractional-Order Sliding Mode Control Method for a Class of Integer-Order Nonlinear Systems
Previous Article in Special Issue
A Novel Decomposition Method for Manufacture Variations and the Sensitivity Analysis on Compressor Blades
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

by
Ernesto Sanchez-Laulhe
1,*,
Ramon Fernandez-Feria
1 and
Anibal Ollero
2
1
Fluid Mechanics Group, University of Malaga, Dr Ortiz Ramos S/N, 29071 Malaga, Spain
2
GRVC Robotics Laboratory, University of Seville, Camino de los Descubrimientos S/N, 41092 Seville, Spain
*
Author to whom correspondence should be addressed.
Aerospace 2022, 9(10), 617; https://doi.org/10.3390/aerospace9100617
Submission received: 14 September 2022 / Revised: 28 September 2022 / Accepted: 14 October 2022 / Published: 18 October 2022
(This article belongs to the Special Issue Aerodynamics Design)

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.

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,8,9]. Such models combine estimations of the aerodynamic forces for high angles of attack, 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 Theodor-sen [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,29,30,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.
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, on-board 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 (scaled with the wing’s chord length) as the 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 2 presents 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 3 shows some numerical results and the experimental validation of the model. Section 4 explains 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 ˜ 0 e i ω t ˜ ,
where ω and h ˜ 0 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].

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:
U c = m g π ρ S δ t , L c = c 2 , t c = 1 ω ,
where U c and t c 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
M k 0 U ˙ = U 2 C T * C D * L i * Λ C D t * δ t sin ( γ ) ,
M k 0 U γ ˙ = U 2 C L * + Λ C L t * δ t cos ( γ ) ,
M k 0 2 θ ¨ = χ M 2 U 2 l w C L * cos ( α ) ( C T * C D * ) sin ( α ) + l t Λ C L t * cos ( α ) + C D t * sin ( α ) h w C L * sin ( α ) + C T * C D * cos ( α ) ,
and the wing’s angle of attack is computed as α = θ γ . Equations (3)–(5) are written in the trajectory frame, with x T in the direction of the velocity. The non-dimensional parameters appearing in Equations (3)–(5) that differ from those used in [27] are
M = 2 m π ρ S c , χ = π ρ S c 2 3 2 I y , k 0 = ω c 2 U c ,
The (modified) force coefficients C L * , C T * , C D * , modeled in Section 2.2 are defined as
C L * = L π ρ U b 2 c , C T * = T π ρ U b 2 c , C D * = D π ρ U b 2 c ,
applied at the wing aerodynamic center. The Lighthill number L i * = S b S C D b * , 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 ( l w , h w ) in the 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 h 0 and the reduced (dimensionless) frequency, k = k 0 U , which is affected by the dynamics through the non-dimensional velocity U in relation to the constantly reduced frequency k 0 , defined in (6).
Since we assume small flapping amplitudes, we use the well-known results from two-dimensional, 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 α 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 Aerospace 09 00617 i001/(Aerospace 09 00617 i001 + 2), where Aerospace 09 00617 i001 = b 2 / S is the aspect ratio of the wing. On the other hand, Smits [16] proposes a discontinuous factor for the added mass terms, with Aerospace 09 00617 i001 / 2 for Aerospace 09 00617 i001 < 2 and just unity for Aerospace 09 00617 i001 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 * = C L α α + C L h e i t + C ¯ L h e i t U h 0 + C L α ˙ U α ˙ + C L U ˙ U 2 α U ˙ + C L θ ˙ U θ ˙ ,
where we see the harmonic function related to the heaving motion multiplying h 0 . The different coefficients C L . . . 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 * = C T c U 2 h 0 2 + C T α α 2 + C T h α e i t + C ¯ T h α e i t U h 0 α + C T h e 2 i t + C ¯ T h e 2 i t U 2 h 0 2 + C T h α ˙ e i t + C ¯ T h α ˙ e i t U 2 h 0 α ˙ + C T θ ˙ α ˙ U 2 θ ˙ α ˙ + C T α ˙ α U α ˙ α + C T h θ ˙ e i t + C ¯ T h θ ˙ e i t U 2 h 0 θ ˙ + C T θ ˙ α U θ ˙ α + C T θ ˙ U 2 θ ˙ 2 + C T U ˙ α U 2 α 2 U ˙ .
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 C T . . . (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 L t * = C L t α α δ t + C L t θ ˙ U θ ˙ + C L t α ˙ U α ˙ = t 4 α δ t t 4 k 0 l t U θ ˙ + 3 t 8 k 0 U α ˙ ,
Finally, drag from all parts of the vehicle has to be considered. Body drag, represented by the modified Lighthill number L i * , 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 L t * , which is computed from finite wing theory. In summary,
C D * = C D 0 * + C D i * , C D i * = 2 C L * 2 ,
C D t * = C D 0 t * + C D i t * , C D i t * = 2 C L t * 2 t .

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 2 are numerically integrated. We use the values of the parameters from the prototype described in [5], which are shown in Table 1 for three flapping frequencies f relevant for bird-scale flight (actually, only k 0 depends on f). Note that they are of order unity, except for the friction drag ( L i * ).
With these parameters, we proceed to numerically integrate the system of differential Equations (3)–(5). Figure 2 shows 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 h 0 = 0.1 . To better compare the three cases, velocities have been normalized with their final values U f , 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.
Figure 3 shows the evolution of the angular variables θ and α with time for a frequency of 5 Hz. 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,29,30,31]. Notice also that the aerodynamic models have no limitation in relation to the frequencies.
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 h 0 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 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).
The prototype has a flapping amplitude h 0 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,29,30,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 C L t θ ˙ 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 5 shows 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 2 s 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 U and θ .

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, ϵ h 0 1 . Before that, some assumptions have to be made regarding the order of magnitude, in terms of ϵ , of the non-dimensional variables and the remaining parameters for a typical ornithopter, which are enumerated next.
  • Small flight angles: γ θ α δ t ϵ . 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 ϵ . Small γ and θ mean that aggressive climbing trajectories are not considered. Moreover, small flapping amplitude implies a limited thrust, which is proportional to ϵ 2 (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 L t * ϵ . 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 * ϵ 2 . Clearly from the Formulation (9), with products of the aerodynamic angles and the flapping amplitude.
  • Drag coefficients: C D * C D t * L i * ϵ 2 . 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: M 2 χ M k 0 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 t c ). 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 as Δ U / τ c 1 and Δ γ / τ c 1 , the order of magnitude of the time derivatives, where τ c 1 is the transient timescale in non-dimensional time t, one has
M k 0 Δ U τ c 1 U ( Δ U ) ϵ 2 + ϵ Δ γ ,
M k 0 U Δ γ τ c 1 U ( Δ U ) ϵ + ϵ γ Δ γ ,
where it has been taken into account that the leading terms with U 2 cancel with the corresponding gravity terms so that the velocity terms going with the aerodynamic coefficients in this transient phase are of the order of U Δ U . Since U 1 by the choice of U c and the characteristic time has to be of the same order in both equations (the numerical solutions plotted in Figure 2 and Figure 3 show a similar oscillatory convergence for U and for the angular variables), it necessarily implies that Δ U is of the same order as Δ γ . Thus,
M k 0 1 τ c 1 ϵ 2 + ϵ ,
M k 0 1 τ c 1 ϵ + ϵ 2 ,
and, comparing the leading terms, τ c 1 = M k 0 / ϵ . These leading terms are related to the gravity force and to the effect of the velocity variations on the lift force. However, the damping terms related to the aerodynamic drag are negligible in this time scale. Thus, we have to consider yet another (slower) characteristic time to account for the drag term in (13),
M k 0 Δ U τ c 2 U ( Δ U ) ϵ 2 ,
obtaining the third timescale τ c 2 = M k 0 / ϵ 2 . Aerodynamic coefficient variations have not been considered as, even though they can affect the transient phase, they will not change the corresponding timescales.
Therefore, the complexity of the present problem needs three timescales to obtain a meaningful perturbation solution: t, τ 1 = t / τ c 1 , and τ 2 = t / τ c 2 . Non-dimensional time derivatives are thus defined as the sum of three terms
d d t = t + ϵ M k 0 τ 1 + ϵ 2 M k 0 τ 2 .
The state variables can then be expressed as an asymptotic expansion in powers of ϵ ,
U = U 0 ( t , τ 1 , τ 2 ) + ϵ U 1 ( t , τ 1 , τ 2 ) + ϵ 2 U 2 ( t , τ 1 , τ 2 ) +
θ = ϵ θ 1 ( t , τ 1 , τ 2 ) + ϵ θ 2 ( t , τ 1 , τ 2 ) +
α = ϵ α 1 ( t , τ 1 , τ 2 ) + ϵ α 2 ( t , τ 1 , τ 2 ) +
Note that α is considered here as a state variable instead of the flight path angle γ to simplify the integration of the aerodynamic forces in the dynamic equations. γ can be expressed as the difference between the pitch angle and the angle of attack:
γ = ϵ θ 1 ( t , τ 1 , τ 2 ) α 1 ( t , τ 1 , τ 2 ) + ϵ 2 [ θ 2 ( t , τ 1 , τ 2 ) α 2 ( t , τ 1 , τ 2 ) ] +
The dependence with each timescale of the different functions in the expansions is derived in Section 4.2 below. Additionally, since the first timescale t is related to harmonic oscillations with the flapping frequency, the different functions in Equations (19)–(21) will be expanded in Fourier series with the necessary number of terms to account for the aerodynamic forces (8) and (9). The harmonic terms are presented in their exponential form. To simplify the expressions, their conjugates are omitted for clarity purposes, but obviously they are included when implemented in the equations. In order to not repeat terms, the Fourier series coefficients of U , θ , α will be named with the capital letters V , T , A , with subscripts following the order in the expansion. For instance, the velocity terms will be expanded as follows,
U 0 = V 0 ( τ 1 , τ 2 ) + V 1 ( τ 1 , τ 2 ) e i t + + V n ( τ 1 , τ 2 ) e i n t
U 1 = V n + 1 ( τ 1 , τ 2 ) + V n + 2 ( τ 1 , τ 2 ) e i t + + V n + k + 1 ( τ 1 , τ 2 ) e i k t
where the number of terms in each order ( n , k , ) is given by the expansion of the equations in the next Section 4.2. These Fourier coefficients are complex, and to also include the phase of the oscillations, the real value needed for the evolution of the system is obtained by taking into account the corresponding complex conjugate.

4.2. Perturbation Solution

Once the expansions (19)–(21) and the time derivatives (18) have been defined, we apply the standard multiple-scales perturbation method to simplify the problem by separating Equations (3)–(5) into a set of simpler equations by equating similar powers of ϵ [38].
In addition to the general perturbation scheme, we must consider the expansion in a Fourier series. Therefore, each equation obtained for a particular power of ϵ can be divided into different equations with the harmonic functions as a common factor, obtaining equations with different harmonic terms ( e 0 , e i t , e 2 i t , …).
This procedure provides a set of simple algebraic equations for each power of ϵ and for each specific harmonic term. At the lowest order O ( ϵ 0 ) = O ( 1 ) , just one term remains, belonging to Equation (3),
U 0 t = 0
meaning that the lowest order of the velocity does not depend on the fast time variable t of the oscillations. This result agrees with the numerical solution in Figure 2, where we saw that the oscillations of the velocity with the flapping frequency have very small amplitude.
Following the mathematical development, at the next order O ( ϵ ) , we obtain the following equations:
U 0 τ 1 + U 1 t = 0
M k 0 U 0 θ 1 t α 1 t = U 0 2 C L α α 1 + Λ C L t α α 1 δ t * + U 0 C L h e i t + C L θ ˙ + Λ C L t θ ˙ θ 1 t + C L α ˙ + Λ C L t α ˙ α 1 t δ t *
M k 0 2 M 2 χ 2 θ 1 t 2 = U 0 2 l w C L α α 1 + Λ l t C L t α α 1 δ t * + U 0 l w C L h e i t + l w C L θ ˙ + Λ l t C L t θ ˙ θ 1 t + l w C L α ˙ + Λ l t C L t α ˙ α 1 t
From the first Equation (26), it follows that U 1 is also independent of the fast time t. Note that, according to (25), terms with U 0 / t disappear from Equations (27) and (28).
Now we separate Equations (26)–(28) by harmonic terms. First, non-oscillatory terms in Equations (26)–(28) generate a system of equations for the coefficients associated with the non-harmonic terms in the Fourier series. From (26), only U 0 / τ 1 = 0 remains, so that the lowest order velocity does not depend on τ 1 . With the other two equations, we have a system for the two variables V 0 and A 0 , without any temporal derivative. Therefore, A 0 and V 0 are, in fact, constants, independent of any time variable, and are given by
A 0 = l t Λ C L t α δ t * l w C L α + l t Λ C L t α
V 0 = δ t * C L α A 0 + C L t α Λ A 0 δ t * .
C L α depends on the velocity V 0 through the reduced frequency in Theodorsen’s function (see Equation (A1)). Note that only the lowest-order velocity has to be considered in Theodorsen’s function to be consistent with the present order of the expansion. In any case, this is a quite simple couple of algebraic equations for A 0 and V 0 . These and all the other algebraic equations that will appear in the development of the perturbation method are solved numerically using Matlab’s function fsolve.
The number of terms in the Fourier series of each variable is limited by the equations, as discussed above. In this order of ϵ , Equation (26) has no oscillating terms for the variable U 1 . Therefore, this variable does not depend on t, and we have that the two lower orders of the velocity only have the non-harmonic term from the Fourier series. Then, Equations (27) and (28) have oscillating terms with just the first mode, which leads to θ 1 and α 1 having only one oscillating term. The corresponding functions are obtained from the following system,
b l u e i M k 0 C L θ ˙ + Λ C L t θ ˙ T 1 i M k 0 + i C L α ˙ + Λ C L t α ˙ + V 0 C L α + Λ C L t α A 1 = C L h
k 0 2 χ + i V 0 l w C L θ ˙ + Λ l t C L t θ ˙ T 1 V 0 2 l w C L α + Λ l t C L t α + i V 0 l w C L α ˙ + Λ l t C L t α ˙ A 1 = l w V 0 C L h
We see that the system of Equations (31) and (32) does not have any derivative with τ 1 and τ 2 , so all terms are constants and, consequently, so are A 1 and T 1 . Note that this is a complex algebraic system, as complex unit i appears with derivatives with variable t. Coefficients A 1 and T 1 are thus complex, giving information about the amplitudes of the oscillations | T 1 | and | A 1 | and their respective phases ϕ T 1 and ϕ A 1 : T 1 = T 1 e i ϕ T 1 and A 1 = A 1 e i ϕ A 1 .
In addition, there is also a transient mode with time scale t from Equations (26)–(28), given by
M k 0 C L α ˙ + Λ C L t α ˙ α 1 t = C L θ ˙ + Λ C L t θ ˙ M k 0 θ 1 t + V 0 C L α + Λ C L t α α 1
k 0 2 χ 2 θ 1 t 2 V 0 l w C L α ˙ + Λ l t C L t α ˙ α 1 t = V 0 l w C L θ ˙ + Λ l t C L t θ ˙ θ 1 t + V 0 2 l w C L α + Λ l t C L t α α 1
From these linear equations, a fast mode can be easily obtained. However, these transients are not of interest. To not alter the dynamics of the ornithopter, we just have to ensure that these fast solutions are stable.
Summarizing, order ϵ of the equation provides information about the lowest order of the three variables. For the velocity, we have the (constant) only term of U 0 , knowing also that U 1 does not have any harmonic oscillation with t. For the pitch angle, we have the harmonic oscillations with t of θ 0 . Finally, all the terms of α 0 are already known.
There is not enough information to obtain the evolution of the system yet, so we have to analyze the order ϵ 2 of the equations. Each Fourier coefficient in this order is expected to have an asymptotic value and a transient term. To simplify the system further, each coefficient is decomposed as X ( τ 1 , τ 2 ) = X s + X t ( τ 1 , τ 2 ) from this point onwards to divide it into two separate systems of equations. Firstly, for the constant terms (with subscripts s), the following analytical expressions are obtained
T 0 s = A 0 + V 0 2 δ t * C T α Λ C D t α A 0 2 + 2 A 1 A ¯ 1 Λ C D t α δ t δ t * A 0 C D 0 * L i * Λ C D 0 t * + V 0 δ t * C T h α A ¯ 1 + C ¯ T h α A 1 + C T θ ˙ α i T 1 A ¯ 1 i T ¯ 1 A 1 + i δ t * C ¯ T h α ˙ A 1 C T h α ˙ A ¯ 1 + C T θ ˙ α ˙ δ t * A 1 T ¯ 1 + A ¯ 1 T 1 + i δ t * C ¯ T h θ ˙ T 1 C T h θ ˙ T ¯ 1 + 2 C T θ ˙ T 1 T ¯ 1 δ t * + C T c δ t *
A 2 s = h w l w C L α + Λ l t C L t α C T α + C L α A 0 2 + 2 A 1 A ¯ 1 C D 0 * + 1 V 0 ( C T h α + C L h ) A ¯ 1 + ( C ¯ T h α + C ¯ L h ) A 1 + ( C T θ ˙ α + C L θ ˙ ) i T 1 A ¯ 1 i T ¯ 1 A 1 + 1 V 0 2 i C ¯ T h α ˙ A 1 C T h α ˙ A ¯ 1 + C T θ ˙ α ˙ A 1 T ¯ 1 + A ¯ 1 T 1 + i C ¯ T h θ ˙ T 1 C T h θ ˙ T ¯ 1 + 2 C T θ ˙ T 1 T ¯ 1 + C T c
V 1 s = A 2 s V 0 2 C L α + Λ C L t α C L α A 0 + C L t α Λ A 0 δ t *
which are easily computed using the terms defined or obtained above. These constant terms provide the asymptotic values to be reached by each state variable after the transient mode at the corresponding order of the expansion.
For the transient terms, which depend on the time variables τ 1 and τ 2 , we obtain
d V 1 t ( τ 1 , τ 2 ) d τ 1 = δ t * T 0 t ( τ 1 , τ 2 )
d T 0 t ( τ 1 , τ 2 ) d τ 1 = 2 V 1 t ( τ 1 , τ 2 ) C L α A 0 + C L t α Λ A 0 δ t * 1 C L θ ˙ + Λ C L t θ M k 0 + l w C L θ ˙ + l t Λ C L t θ ˙ M k 0 C L α + Λ C L t α l w C L α + l t Λ C L t α
A 2 t ( τ 1 , τ 2 ) = d T 0 t ( τ 1 , τ 2 ) d τ 1 l w C L θ ˙ + Λ l t C L t θ ˙ M k 0 V 0 l w C L α + Λ l t C L t α
which is a system of two linear differential equations and a third algebraic expression for the evolution of the angle of attack. The solution to this system of ordinary differential equations has the form X ( τ 1 , τ 2 ) = e λ τ 1 X ( τ 2 ) , where λ turns out to be imaginary, λ = i ω 1 , representing a purely oscillating transient phase with frequency
ω 1 = 2 δ t * C L α A 0 + C L t α Λ A 0 δ t * 1 C L θ ˙ + Λ C L t θ M k 0 + l w C L θ ˙ + l t Λ C L t θ ˙ M k 0 C L α + Λ C L t α l w C L α + l t Λ C L t α
related to the mean lift of the ornithopter ( C L α ) and the aerodynamic coefficient associated with the pitch angular velocity. It is observed that the frequency of this transient phase is not affected by the flapping frequency, in agreement with the numerical results shown in Section 3.1.
Given the form of the solution, the functions of τ 2 that remain to be solved are complex, as they also include the information related to the phase of the oscillations. However, we can write them as real variables computing their phases separately and by using Equation (38) and the initial values U i and T i to yield
i ω 1 e i ϕ V V 1 t ( τ 2 ) = e i ϕ T δ t * T 0 t ( τ 2 )
ϕ V π 2 = ϕ T , V 1 t ( τ 2 ) = δ t * T 0 t ( τ 2 ) ω 1
tan ϕ T = ω 1 U i δ t * T i
Equation (43) gives the relations between phases and amplitudes of the two variables obtained from Equation (42). The phase is determined by the initial conditions, as in Equation (44), and the phase shift is always 90 deg. This transient mode is common for aircraft and is known as phugoid mode. For the evolution with τ 2 , just another equation is needed, as the relation between V 1 t ( τ 2 ) and T 0 t ( τ 2 ) is already known.
Regarding the harmonic terms at the order ϵ 2 of Equations (3)–(5), both first and second-frequency terms, i.e., proportional to e i t and e 2 i t , appear in the equations. Hence, the three variables U 2 , θ 1 and α 1 have three terms in the Fourier series expansion. Considering Equation (3), the two resulting velocity coefficients are
V 3 = i M k 0 V 0 2 C T α Λ C D t α 2 A 0 A 1 Λ C D t α δ t δ t * A 1 + V 0 A 0 C T h α + i C T α ˙ α A 1 + i C T θ ˙ α T 1 δ t * T 1 + δ t * A 1
V 4 = i 2 M k 0 V 0 2 A 1 2 C T α Λ C D t α + V 0 A 1 C T h α + i C T α ˙ α A 1 + i C T θ ˙ α T 1 + C T h + i C T h α ˙ A 1 C T θ ˙ α ˙ A 1 T 1 + i C T h θ ˙ T 1 C T θ ˙ T 1 2
Note that these lowest-order oscillations do not depend on the transient timescale, and they are constant during the entire evolution of the system. Furthermore, unlike the lowest-order oscillations for the pitch angle and the angle of attack considered above, both first and second-frequency oscillations appear in the velocity.
The oscillatory terms from (4) and (5) lead to the coefficients A 3 and T 3 for the harmonic term e i t , and A 4 and T 4 for e 2 i t . However, the analysis is algebraically more involved, as summarized in Appendix B. As it happens with coefficients V 1 , T 0 and A 2 , A 3 and T 3 contain steady and transient terms.
Finally, to obtain T 0 t ( τ 2 ) and V 1 t ( τ 2 ) , order ϵ 3 of Equations (3)–(5) has to be considered. In this case, variables of higher order with the frequency ω 1 appear. When combining the three equations for the transient non-harmonic terms, a single differential equation is obtained for the evolution of T 0 t ( τ 2 ) (or alternatively V 1 t ( τ 2 ) ). The result is written as T 0 t ( τ 2 ) = T i e ( i ω 2 ξ ) τ 2 . The expressions for ω 2 and ξ are found in Appendix B.

4.3. Comparison between Analytical and Numerical Results

The theoretical approach developed allows a complete characterization of flapping forward flight. In contrast to other previous works [18,39], it also takes into account oscillations in state variables due to flapping wing movement, with a simplified model appropriate for on-board computers. This section presents a comparison between analytical and numerical results, showing their excellent agreement.
In summary, the analytical results obtained in the previous section are the following:
  • For the large-time asymptotic values of the velocity and the pitch angle, we consider terms up to order ϵ 2 , while for the angle of attack, the approximation is even better, up to order ϵ 3 .
  • Transient terms are computed up to order ϵ for the velocity and the pitch angle, and to order ϵ 2 for the angle of attack.
  • Harmonic oscillations with the flapping frequencies are considered up to order ϵ 2 for the three variables, including asymptotic and transient values.
Mathematically, the analytical solution used in the comparison is (the different constants and functions are obtained in Section 4.2 or Appendix B):
U V 0 + ϵ V 1 s + V 1 t ( τ 1 , τ 2 ) + ϵ 2 V 2 s + V 3 e i t + V 4 e 2 i t ,
θ ϵ T 0 s + T 0 t ( τ 1 , τ 2 ) + ϵ T 1 e i t + ϵ 2 T 2 s + T 3 ( τ 1 , τ 2 ) e i t + T 4 ( τ 1 , τ 2 ) e 2 i t ,
α ϵ A 0 + A 1 e i t + ϵ 2 A 2 s + A 2 t ( τ 1 , τ 2 ) + A 3 ( τ 1 , τ 2 ) e i t + A 4 ( τ 1 , τ 2 ) e 2 i t + ϵ 3 A 5 s .
Physically, the large-time permanent state consists of an oscillation of the variables in the fast timescale t (i.e., with the flapping frequency and some of its higher harmonics) around a steady value. The perturbation approach gives a low amplitude of the airspeed oscillations compared to its average value. The transient phase occurs in the slowest timescale τ 2 , while the period of the oscillations in the convergence towards the permanent state is of the order of the intermediate timescale τ 1 .
We present three different comparisons between the analytical perturbation solution and the numerical results: for the permanent state, for the transient phase and for the resulting trajectories. All the reported computations are for a frequency of 5 Hz, an amplitude ϵ = h 0 = 0.1 , and a tail deflection of 4 . Similar results have been obtained for frequencies between 2 and 7 Hz, amplitudes between 0.05 and 0.3 , and tail deflections between 0 . 5 and 10 . Limits are fixed by aerodynamic theory for both amplitude and tail deflection. The frequency range is based on existing prototypes [5].
Figure 6 shows a good agreement between the numerical results and the analytical solution (47)–(49) for the large-time permanent state. Differences are of one order higher than that of the last term considered in (47)–(49), so we can conclude the validity of the analytical method for the given conditions. Figure 6a also assesses some of the hypotheses proposed for the simplified method, namely the order unity of the asymptotic airspeed and the order ϵ 2 of its oscillation amplitude. The minimal error of the perturbation solution for U does not suppose a significant difference in the trajectory followed by the UAV, and for long trajectories, it does not cause a larger error than that originated by inaccuracies in the measurement of the UAV characteristics. On the other hand, Figure 6b corroborates that the time-averaged value of the pitch angle is of the same order as the amplitude of its oscillations (their actual values in this case are even smaller than O ( ϵ ) ), showing an excellent agreement between the analytical prediction and the numerical solution.
Figure 7 compares the transient phase for both the velocity and the pitch angle. Note the similarity in the shape of both transient phases, indicating the good adjustment of the timescales. However, the agreement with the numerical results of the transient phase is not as good as the one observed for the permanent flapping state in Figure 6. We also see how the orders considered in Section 3.1 are correct, with the transient velocity one order lower than its asymptotic value. On the other hand, the transient component of the pitch angle is of order ϵ , even though its asymptotic value is smaller in this particular case. Notice also that the amplitudes of the pitch oscillations remain in the adequate range for the reasonable application of the aerodynamic models.
Figure 8 shows 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.

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 real-time 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 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,29,30,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.

Appendix A. Aerodynamic Coefficients in Equations (9) and (10)

C L α = + 2 F k ,
C L h = k 0 2 k 2 + + 2 G ( k ) F ( k ) i ,
C L α ˙ = C L U ˙ = k 0 2 ,
C L θ ˙ = + 2 F ( k ) k 0 l w 1 ,
where F ( k ) and G ( k ) are the real and imaginary parts of Theodorsen’s function [10].
C T c = k 0 2 π G 1 k + 2 ,
C T α = 2 π G 1 k + 2 C L α ,
C T h α = k 0 π F 1 k + 2 i G 1 k + 2 C L h ,
C T h = i k 0 2 2 π F 1 k + i G 1 k + 2 ,
C T h α ˙ = i k 0 2 4 k 0 2 2 + 4 i F 2 k + G 2 k ,
C T θ ˙ α ˙ = k 0 2 2 l w 1 2 k 0 2 l w 1 F 2 k + 2 ,
C T α ˙ α = k 0 2 + k 0 F 2 k + 2 C L α ˙ ,
C T h θ ˙ = k 0 2 π i G 1 k + 2 l w 1 + k 0 2 π i G 1 k F 1 k + 2 l w 1 2 ,
C T θ ˙ α = 2 k 0 π G 1 k + 2 l w 1 + 2 k 0 π G 1 k + 2 l w 1 2 C L θ ˙ ,
C T θ ˙ = 2 k 0 2 π G 1 k + 2 l w 1 l w 1 2
C T U ˙ α = C L U ˙ ,
where F 1 ( k ) and G 1 ( k ) are the real and imaginary parts of the modified function defined in [14], while F 2 ( k ) and G 2 ( k ) are given by
F 2 ( k ) = G ( k ) + 2 π k G 1 ( k ) F 1 ( k ) , G 2 ( k ) = F ( k ) + 2 π k F 1 ( k ) + G 1 ( k ) .

Appendix B. Higher Order Perturbation Terms

System of algebraic equations for T 3 and A 3 :
M 1 T 3 A 3 = b 1
being
M 1 ( 1 , 1 ) = i M k 0 C L θ ˙ Λ C L t θ ˙
M 1 ( 1 , 2 ) = i M k 0 + C L α ˙ + Λ C L t α ˙ V 0 C L α + Λ C L t α
M 1 ( 2 , 1 ) = k 0 2 χ + i V 0 l w C L θ ˙ + l t Λ C L t θ ˙
M 1 ( 2 , 2 ) = V 0 2 l w C L α + l t Λ C L t α + i V 0 l w C L α ˙ + l t Λ C L t α ˙ ,
and
b 1 ( 1 ) = 2 V 1 C L α + Λ C L t α A 1 + V 1 V 0 C L h i T 1 M k 0 C L θ ˙ Λ C L t θ ˙ + A 1 M k 0 + C L α ˙ + Λ C L t α ˙
b 1 ( 2 ) = 2 V 0 V 1 A 1 l w C L α + l t Λ C L t α + V 1 i T 1 l w C L θ ˙ + l t Λ C L t θ ˙ + i A 1 l w C L α ˙ + l t Λ C L t α ˙ + l w V 1 C L h V 0 h w A 0 V 0 C T α + C L α 2 A 1 + C T h α + C L h + C L α ˙ + C T α ˙ α i A 1 + C L θ ˙ + C T θ ˙ α i T 1
Note that T 3 and A 3 have a transient term. Its solution can be written as
T 3 = T 3 s + V 1 t ( τ 1 , τ 2 ) z T 3
A 3 = A 3 s + V 1 t ( τ 1 , τ 2 ) z A 3
where z T 3 and z A 3 are complex and constant.
System of algebraic equations for T 4 and A 4 :
M 2 T 4 A 4 = b 2
being
M 2 ( 1 , 1 ) = 2 i M k 0 C L θ ˙ Λ C L t θ ˙
M 2 ( 1 , 2 ) = 2 i M k 0 + C L α ˙ + Λ C L t α ˙ V 0 C L α + Λ C L t α
M 2 ( 2 , 1 ) = 4 k 0 2 χ + 2 i V 0 l w C L θ ˙ + l t Λ C L t θ ˙
M 2 ( 2 , 2 ) = V 0 2 l w C L α + l t Λ C L t α + 2 i V 0 l w C L α ˙ + l t Λ C L t α ˙ ,
and
b 2 ( 1 ) = 0
b 2 ( 2 ) = h w V 0 2 C T α + C L α + V 0 A 1 C T h α + C L h + C L α ˙ + C T α ˙ α i A 1 + C L θ ˙ + C T θ ˙ α i T 1 + C T h + i C T h α ˙ A 1 C T θ ˙ α ˙ A 1 T 1 + i C T h θ ˙ T 1 C T θ ˙
Expressions for T 2 s , A 5 s and V 2 s :
T 2 s = A 2 s + 2 V 0 V 1 δ t * C T α Λ C D t α A 0 2 + 2 A 1 A ¯ 1 Λ C D t α δ t δ t * A 0 C D 0 * L i * Λ C D 0 t * + V 1 s δ t * C T h α A ¯ 1 + C ¯ T h α A 1 + i C T θ ˙ α T 1 A ¯ 1 T ¯ 1 A 1 + V 0 2 δ t * C T α Λ C D t α A 0 A 2 s + A 1 A ¯ 3 s + A ¯ 1 A 3 s Λ C D t α δ t δ t * A 2 s + V 0 δ t * C T h α A ¯ 3 s + C ¯ T h α A 3 s + i C T θ ˙ α T 3 s A ¯ 1 T ¯ 3 s A 1 + T 1 A ¯ 3 s T ¯ 1 A 3 s + i δ t * C ¯ T h α ˙ A 3 s C T h α ˙ A ¯ 3 s + C T θ ˙ α ˙ δ t * T 3 s A ¯ 1 T ¯ 3 s A 1 + T 1 A ¯ 3 s T ¯ 1 A 3 s + i δ t * C ¯ T h θ ˙ T 3 s C T h θ ˙ T ¯ 3 s + 2 C T θ δ t * T 1 T ¯ 3 s + T ¯ 1 T 3 s
A 5 s = 1 V 0 2 l w C L α + l t Λ C L t α 2 V 0 V 3 l w C L α + l t Λ C L t α A ¯ 1 + 2 V 0 V ¯ 3 l w C L α + l t Λ C L t α A 1 + V 3 l w C ¯ L h i A ¯ 1 l w C L α ˙ + l t Λ C L t α ˙ i T ¯ 1 l w C L θ ˙ + l t Λ C L t θ ˙ + V ¯ 3 l w C L h i A 1 l w C L α ˙ + l t Λ C L t α ˙ i T 1 l w C L θ ˙ + l t Λ C L t θ ˙ + 2 V 0 V 1 s l w C L α + l t Λ C L t α A 2 s h w C T α + C L α A 0 2 + 2 A 1 A ¯ 1 Λ C D t α δ t δ t * A 0 C D 0 * V 1 s h w C L h + C T h α A ¯ 1 + C ¯ L h + C ¯ T h α A 1 + i C L θ ˙ + C T θ ˙ α T 1 A ¯ 1 T ¯ 1 A 1 h w i C ¯ T h α ˙ A 3 s C T h α ˙ A ¯ 3 s + C T θ ˙ α ˙ A 3 s T ¯ 1 + A ¯ 3 s T 1 + A 1 T ¯ 3 s + A ¯ 1 T 3 s + i C ¯ T h θ ˙ T 3 s C T h θ ˙ T ¯ 3 s h w C T θ ˙ T 1 T ¯ 3 s + T ¯ 1 T 3 s + i l w C L U ˙ A ¯ 1 V 3 + A 1 V ¯ 3 + V 0 2 l w C L α 2 + C T α A 0 3 + 6 A 0 A 1 A ¯ 1 + A 0 C D 0 * + V 0 2 l t Λ C L t α 2 C D t α A 0 3 + 6 A 0 A 1 A ¯ 1 + C L t α 2 C D t α δ t A 0 2 + A 1 A ¯ 1 δ t * + A 0 C D 0 * 2 h w C T α + C L α A 0 A 2 s + A 1 A ¯ 3 s + A ¯ 1 A 3 s
V 2 s = V 1 s V 0 + 1 2 V 0 C L α A 0 + Λ C L t α A 0 δ t * i M k 0 V ¯ 3 T 1 A 1 V 3 T ¯ 1 A ¯ 1 V 3 C ¯ L h i A ¯ 1 C L α ˙ + Λ C L t α ˙ i T ¯ 1 C L θ ˙ + Λ C L t θ ˙ 2 V 0 V 3 C L α + Λ C L t α A ¯ 1 V ¯ 3 C L h i A 1 C L α ˙ + Λ C L t α ˙ i T 1 C L θ ˙ + Λ C L t θ ˙ 2 V 0 V ¯ 3 C L α + Λ C L t α A 1 i C L U ˙ A ¯ 1 V 3 A 1 V ¯ 3 2 V 0 V 1 s A 2 s C L α + Λ C L t α V 0 2 A 5 s C L α + Λ C L t α δ t * A 0 2 + A 1 A ¯ 1 + T 0 s 2 + T 1 T ¯ 1 A 0 T 0 s A 1 T ¯ 1 A ¯ 1 T 1
Damping of the transient phase ξ
ξ = C L U ˙ A 0 2 Γ M k 0 1 + Δ l w ω 1 2 Δ 2 Γ M 2 χ C T h α A ¯ 1 + C ¯ T h α A 1 + i C T θ ˙ α T 1 A ¯ 1 T ¯ 1 A 1 + i C ¯ T h α ˙ z A 3 C T h α ˙ z ¯ A 3 2 C T θ ˙ α ˙ z A 3 T ¯ 1 + z ¯ A 3 T 1 + z T 3 A ¯ 1 + z ¯ T 3 A 1 + i C ¯ T h θ ˙ z T 3 + C T h θ ˙ z ¯ T 3 + 2 C T θ ˙ T 1 z ¯ T 3 + T ¯ 1 z T 3 2 V 0 2 C T h α z ¯ A 3 + C ¯ T h α z A 3 i C T θ ˙ α z A 3 T ¯ 1 + z ¯ A 3 T 1 A 1 z ¯ T 3 + A ¯ 1 z T 3 δ t * T 0 s A 0 2 Γ V 0 C T α Λ C D t α A 0 2 + 2 A 1 A ¯ 1 Λ C D t α δ t δ t * A 0 C D 0 * L i * Λ C D 0 t * V 0 C T θ ˙ α A 0 ω 1 2 2 M k 0 V 0 2 2 C T α Λ C D t α i ω 1 A 0 z A 2 + A 1 z ¯ A 3 + A ¯ 1 A 3 s + i ω 1 Λ C D t α δ t δ t * z A 2 + δ t * i ω 1 z A 2 2
and correction of the transient frequency with τ 2 , ω 2
ω 2 = V 1 s Γ ω 1 C L α A 0 + Λ C L t α A 0 δ t * + C L α 1 + Δ l w + Λ C L t α 1 + Δ l t Γ V 0 ω 1 A 2 s V 0 V 1 i z A 4 + C L θ ˙ 1 + Δ l w + Λ C L t θ ˙ 1 + Δ l t V 1 s ω 1 2 Γ M k 0 Δ V 0 Γ ω 1 C T α + C L α A 0 2 + 2 A 1 A ¯ 1 C D 0 * Δ h w 2 Γ ω 1 C L h + C T h α A ¯ 1 + C ¯ L h + C ¯ T h α A 1 + i C L θ ˙ + C T θ ˙ α T 1 A ¯ 1 T ¯ 1 A 1 + i C ¯ T h α ˙ z A 3 C T h α ˙ z ¯ A 3 Δ h w 2 Γ ω 1 C T θ ˙ α ˙ z A 3 T ¯ 1 + z ¯ A 3 T 1 + z T 3 A ¯ 1 + z ¯ T 3 A 1 + i C ¯ T h θ ˙ z T 3 + C T h θ ˙ z ¯ T 3 + 2 C T θ ˙ T 1 z ¯ T 3 + T ¯ 1 z T 3 Δ V 0 2 h w Γ ω 1 C T α + C L α A 0 i ω 1 z A 2 + A 1 z ¯ A 3 + A ¯ 1 z A 3 Δ V 0 h w ω 1 2 Γ M k 0 C T θ ˙ α + C L θ ˙ A 1 V 1 s ω 1 2 Γ Δ V 0 h w 2 Γ ω 1 C T h α + C L h z ¯ A 3 + C ¯ T h α + C ¯ L h z A 3 i C T θ ˙ α + C L θ ˙ z A 3 T ¯ 1 + z ¯ A 3 T 1 A 1 z ¯ T 3 + A ¯ 1 z T 3
with A 2 ( τ 1 , τ 2 ) defined as
A 2 t ( τ 1 , τ 2 ) = z A 2 T 0 t ( τ 2 ) e Φ T + i ω 1 τ 1 ; z A 2 = i ω 1 M k 0 V 0 l w C L θ ˙ + Λ l t C L t θ ˙ l w C L α + Λ l t C L t α
and the constant terms
Γ = V 0 1 C L θ ˙ + Λ C L t θ ˙ M k 0 + l w C L θ ˙ + Λ l t C L t θ ˙ M k 0 Δ
Δ = C L α + Λ C L t α l w C L α + Λ l t C L t α

References

  1. Mueller, T.J. Fixed and Flapping Wing Aerodynamics for Micro Air Vehicle Applications; AIAA: Notre Dame, IN, USA, 2001. [Google Scholar]
  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. [Google Scholar] [CrossRef]
  3. Srigrarom, S.; Chan, W.-L. Ornithopter Type Flapping Wings for Autonomous Micro Air Vehicles. Aerospace 2015, 2, 235–278. [Google Scholar] [CrossRef] [Green Version]
  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. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar] [CrossRef]
  8. Wang, Z.J. Dissecting insect flight. Annu. Rev. Fluid Mech. 2005, 37, 183–210. [Google Scholar] [CrossRef] [Green Version]
  9. Hefler, C.; Kang, C.; Qiu, H.; Shyy, W. Distinct Aerodynamics of Insect-Scale Flight; Cambridge University Press: Cambridge, UK, 2021. [Google Scholar] [CrossRef]
  10. Theodorsen, T. General Theory of Aerodynamic Instability and the Mechanism of Flutter; NACA Technical Report TR 496; NACA: Washington, DC, USA, 1935. [Google Scholar]
  11. Garrick, I.E. Propulsion of a Flapping and Oscillating Airfoil; NACA Technical Report TR 567; NACA: Washington, DC, USA, 1936. [Google Scholar]
  12. DeLaurier, J.D. An aerodynamic model for flapping-wing flight. Aeronaut. J. 1993, 97, 125–130. [Google Scholar]
  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: Bellingham, WA, USA, 2008; Volume 6928, pp. 477–484. [Google Scholar]
  14. Fernandez-Feria, R. Linearized propulsion theory of flapping airfoils revisited. Phys. Rev. Fluids 2016, 1, 084502. [Google Scholar] [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. [Google Scholar] [CrossRef]
  16. Smits, A.J. Undulatory and oscillatory swimming. J. Fluid Mech. 2019, 874, 1–44. [Google Scholar] [CrossRef] [Green Version]
  17. Dietl, J.M.; Garcia, E. Stability in ornithopter longitudinal flight dynamics. J. Guid. Control. Dyn. 2008, 31, 1157–1163. [Google Scholar] [CrossRef]
  18. Dietl, J.M.; Garcia, E. Ornithopter optimal trajectory control. Aerosp. Sci. Technol. 2013, 26, 192–199. [Google Scholar] [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. [Google Scholar] [CrossRef] [Green Version]
  20. Taylor, G.K.; Żbikowski, R. Nonlinear time-periodic models of the longitudinal flight dynamics of desert locusts Schistocerca gregaria. J. R. Soc. Interface 2005, 2, 197–221. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  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. [Google Scholar] [CrossRef]
  22. Grauer, J.A.; Hubbard, J.E., Jr. Multibody model of an ornithopter. J. Guid. Control. Dyn. 2009, 32, 1675–1679. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar] [CrossRef] [Green Version]
  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. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar] [CrossRef] [Green Version]
  30. Fernandez-Feria, R. Note on optimum propulsion of heaving and pitching airfoils from linear potential theory. J. Fluid Mech. 2017, 826, 781–796. [Google Scholar] [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. [Google Scholar] [CrossRef]
  32. Sánchez-Rodríguez, J.; Raufaste, C.; Argentina, M. A minimal model of self propelled locomotion. J. Fluids Struct. 2020, 97, 103071. [Google Scholar] [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. [Google Scholar]
  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. [Google Scholar] [CrossRef] [Green Version]
  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, USA, 1946. [Google Scholar]
  36. Thomas, A.L.R. On the aerodynamics of birds’ tails. Philos. Trans. R. Soc. Lond. Ser. B Biol. Sci. 1993, 340, 361–380. [Google Scholar] [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. [Google Scholar] [CrossRef]
  38. Kevorkian, J.; Cole, J.D. Perturbation Methods in Applied Mathematics; Springer: New York, NY, USA, 1981. [Google Scholar]
  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. [Google Scholar] [CrossRef]
Figure 1. Scheme of forces, axis, and control and state variables. All vectors are on the plane ( X , Z ) of the longitudinal flight.
Figure 1. Scheme of forces, axis, and control and state variables. All vectors are on the plane ( X , Z ) of the longitudinal flight.
Aerospace 09 00617 g001
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 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.
Aerospace 09 00617 g002
Figure 3. Evolution of θ (a) and α (b) with non−dimensional time. Frequency 5 Hz.
Figure 3. Evolution of θ (a) and α (b) with non−dimensional time. Frequency 5 Hz.
Aerospace 09 00617 g003
Figure 4. Ornithopter prototype used in the experiments.
Figure 4. Ornithopter prototype used in the experiments.
Aerospace 09 00617 g004
Figure 5. Comparison between numerical and experimental evolutions. (a) Airspeed; (b) pitch angle.
Figure 5. Comparison between numerical and experimental evolutions. (a) Airspeed; (b) pitch angle.
Aerospace 09 00617 g005
Figure 6. Comparison between theoretical and numerical final state at 5Hz. (a) Airspeed; (b) pitch angle.
Figure 6. Comparison between theoretical and numerical final state at 5Hz. (a) Airspeed; (b) pitch angle.
Aerospace 09 00617 g006
Figure 7. Comparison between theoretical and numerical transient phase at 5Hz. (a) Airspeed evolution; (b) pitch angle evolution.
Figure 7. Comparison between theoretical and numerical transient phase at 5Hz. (a) Airspeed evolution; (b) pitch angle evolution.
Aerospace 09 00617 g007
Figure 8. Comparison between the trajectories at 5 Hz.
Figure 8. Comparison between the trajectories at 5 Hz.
Aerospace 09 00617 g008
Table 1. Non-dimensional flight parameters for three flapping frequencies.
Table 1. Non-dimensional flight parameters for three flapping frequencies.
f M k 0 M M 2 χ Λ l w l t Λ h w Li * Aerospace 09 00617 i001
2 Hz 0.79 2.54 2.12 0.25 0.55 1.16 0.38 0.0048 5.14
5 Hz 1.98 2.54 2.12 0.25 0.55 1.16 0.38 0.0048 5.14
7 Hz 2.77 2.54 2.12 0.25 0.55 1.16 0.38 0.0048 5.14
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

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

AMA Style

Sanchez-Laulhe E, Fernandez-Feria R, Ollero A. Simplified Model for Forward-Flight Transitions of a Bio-Inspired Unmanned Aerial Vehicle. Aerospace. 2022; 9(10):617. https://doi.org/10.3390/aerospace9100617

Chicago/Turabian Style

Sanchez-Laulhe, Ernesto, Ramon Fernandez-Feria, and Anibal Ollero. 2022. "Simplified Model for Forward-Flight Transitions of a Bio-Inspired Unmanned Aerial Vehicle" Aerospace 9, no. 10: 617. https://doi.org/10.3390/aerospace9100617

APA Style

Sanchez-Laulhe, E., Fernandez-Feria, R., & Ollero, A. (2022). Simplified Model for Forward-Flight Transitions of a Bio-Inspired Unmanned Aerial Vehicle. Aerospace, 9(10), 617. https://doi.org/10.3390/aerospace9100617

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop