Aerodynamic Reduced-order Volterra Model of an Ornithopter under High-amplitude Flapping.

The unsteady aerodynamics of ﬂapping low-aspect-ratio ellipsoidal-wings in ornithopters is analyzed and modeled by the use of three dimensional Computational Fluid Dynamics (CFD) simulations. The range of interest is large amplitude, moderate frequency ﬂapping, and low to moderate angles of attack at Reynolds around 10 5 , where autonomous ornithopters like GRIFFIN 1 , are able to perform complex maneuvers such as perching. The results obtained show that the Leading Edge Vortex is produced above a certain Strouhal and angle of attack at downstroke. These aerodynamic loads are compared with the classical analytical models by the frequency response, observing that analytical models based on abscence of viscosity and small perturbations are not appropriate for the range of interest as the hypotheses are not fulﬁlled. Through the 3D CFD aerodynamic loads database, a ﬁnite memory Volterra model is identiﬁed in order to predict the characteristics of forces and moments produced by the ﬂapping wing. This reduced order model depends on the eﬀective angle of attack of the surrogate airfoil located at 70 % of the semi-span at three-quarters chord on the airfoil. This state has been found appropriate for being the one with the greatest regression, comparing the 3D CFD simulations with others that have been carried out in 2D, in agreement with the literature. Finally, a methodology to validate the identiﬁed model without the need of wind tunnel is proposed and validated for lift force. By the use of the aerodynamic forces extracted from ﬂight data, measured by a high accuracy Motion Capture System at diﬀerent ﬂapping wing kinematics, it is concluded that the model provides better estimates than classical analytical models. The structure of the model and its predictability make it possible to use it in control tuning, in addition to being able to append nonlinear or aeroelastic terms using a similar method, also because of the execution time, provide a potential solution for online forces prediction. Bio-inspired UAV, Volterra Model,


Introduction
The efficiency and maneuverability, both in gliding and flapping [1], of bird flight has fascinated scientists throughout history, a fact that can be shown in the diversity of ornithopter prototypes developed, [2], [3], [4]. Emphasis is currently placed on this thanks to several projects in which new capabilities of bio-inspired UAVs are investigated [5]. The optimization of the aerodynamic loads such as lift can enhance the payload of the ornithopters such the GRIFFIN's prototype, figure 1, being able to improve the sensing and computational capabilities and allow the accomplishment of fully autonomous flight and perching. Besides, aerodynamic models are needed to improve designs and for model based control and planning. However aerodynamic modelling of high amplitude flapping wing is a challenging task in terms of accuracy and complexity since involves unsteady viscous phenomena on three dimensional flow. While simple analytical models subject to various hypotheses are practical in the preliminary stages of ornithopter design, they become overly complex when trying to model flow more accurately even more where their hypothesis are not fulfilled. In addition, the physical meaning of the governing variables tends to vanish, making difficult to validate them through experiments. Therefore, CFDbased models provides a useful tool in the middle stages of ornithopter development. Post-processing techniques provide a broad view of the phenomena that occur in the fluid in the UAV's operating range, facilitating the correct choice of the model to identify. Additionally, the validation process becomes more efficient by focusing on the ornithopter's specific operating range.

Classical analytical models and corrections
Unsteady aerodynamic models reflect the dependence of aerodynamic forces with the temporal history of the airfoil states, beyond quasi-steady models that assume that the forces produced depend only on the position of the wing respect to the airflow at each instant of time. Researchers have demonstrated that this aerodynamics do not capture the real phenomenon of flapping wing at intermediate-range Reynolds [6], [7]. Early models developed for unsteady aerodynamics were those of Wagner [8] and Theodorsen [9], for the lift and moment forces in airfoils, derived from the well known linearized potential theory and therefore assuming their hypotheses. The first represents the lift evolution after a change in the effective angle of attack at 3/4 chord by the impulse response and the convolution integral. Since the theory was linearized, the contributions of the various aerodynamic phenomena, added mass and circulation, could be calculated separately and added later. On the other hand, Theodorsen derived a model for the frequency response of the same problem, and highlighted the relevance of the reduced frequency in performance. Both models were later shown to be equivalent [10].
A year later, Garrick [11] made use of the linearized potential theory and the Theodorsen function to derive a compact model of the propulsive forces caused by an airfoil performing small oscillatory heaving and pitching. Furthermore this theory has been completed recently by Feria [12], with the use of the Vortical Impulse Theory within the potential Flow framework a new complex function is proposed.
Assuming that a flapping wing's airfoil in heaving could be interpreted as an infinite wing in such movement, surrogate models found similarities between the performance of the airfoils at given points along the span and that of the three-dimensional flapping wing. Based on Wagner's Indicial Response, Jones [13] proposed a correction of the Theodorsen's Theory with the aspect ratio for threedimensional wings. The results shown that the step response of finite ellipsoidal and infinite wings are equivalent at the beginning, and differ as the wake develops, reaching much lower lift values for wings of small AR. Subsequently, some corrected models based on the Lifting Line Theory were derived from the Wager's step response in order to model arbitrary displacements of the wing. The method consists on the discretization of a three-dimensional wing through two-dimensional airfoils, assuming that the hypotheses of the classical models are fulfilled and adding the contributions of each airfoil separately. Delaurier [14] introduced an efficiency coefficient of the leading edge suction proposed by Garrick, due to the viscous effects, as well as a prediction of the static stall of the wing airfoils.

State of the art: Low Reynolds and high amplitude modelling
At present, various investigations have been carried out on the aerodynamics of flapping wings, as well as on aeroelasticity, which is being studied as a potentially performance-enhancing phenomenon of flapping wing systems at low Reynolds. Moreover, there is a marked tendency in flap wing research to study insect and bird flights separately. While insects generally operate at very low speeds through high beat frequencies and complex wing kinematics, where viscous phenomena mainly occur [15], [16], [17], birds fly at higher speeds, where circulatory effects are not negligible [18]. A review of models proposed for insects till the date can be found at [19].
Although classical models based on the lifting line and potential theory are still widely used for their simplicity and physical meaning [20], due to the development of Computational Fluid Dynamics and the accessibility to new experimental techniques, the underlying phenomena at intermediate-range Reynolds has been thoroughly studied, both in airfoils and three dimensional wings. Several researchers agree that the Leading Edge Vortex (LEV) is the cause of hyper-lift effects in flapping wings at moderate Strouhal numbers in the Reynolds range of 10 4 − 10 5 , that contemplates insects and vertebrates [21]. Smoke visualization experiments and CFD calculations agree the statement [22], [6], [23]. Reference [24] presents a recent review on the state of the art of the aerodynamics of flapping airfoils, emphasizing the parameters that govern performance and revealing the most relevant conclusions of the numerical studies and experiments that have been carried out. One of its conclusions is that a further study of the performance of 3D flapping wing aerodynamics is necessary to better understand the flight of birds in order to achieve efficient real applications.
Due to the complexity increase, Reduced Order Models adapted to the operating range of interest have been identified, base on numerical and experimental data, and mainly from a dynamical system perspective in recent years. Brunton et. al [25] proposed a semi-empirical state-space representation of Theodorsen's lift model, by identifying the circulatory and added mass parameters empirically at intermediate-range Reynolds flows. Zakaria et. al. [26] [27] identified by wind tunnel experiments the empirical frequency response of a plunging wing at high angles of attack, resulting in a lift enhancement compared to Theodorsen's prediction at certain frequency. Also identified a frequency response model of the lift at high angles of attack, by linearization at each angle. However, the model is identified by minimizing the magnitude error, but the phase is not optimized. Boutet et. al. [28] extended the Wagner theory to three dimensional wings by using the lifting line theory, maintaining the small perturbations hypothesis. Taha et. al [29] included to Wagner's Indicial Respose the effects of the stall modeled by the lift-alpha slope, using the Duhamel's principle to obtain a reduced order model for insects at high reduced frequencies. Currently [30], the author has shown and analyzed the differences of the frequency response between the Theodorsen and Wagner's model and the extracted from laminar Navier Stokes simulation of a airfoil. The main conclusions are that, at high frequencies, the frequency response slope of a flapping airfoil differs from the classical theories, being different that zero, and proposed this effect as a "viscous circulation". His conclusions are supported by showing that the classification between circulatory and non-circulatory lift is arbitrary, and also highlights that due to the lift circulation term of Theodorsen differs from vorciticy lift, the Theodorsen transfer function must contain some added mass effect since there are instantaneous lift in his formulation.
Moreover, Silva et. al. [31], [32] proposed the Volterra theory to high velocity unsteady aerodynamic responses analysis. Cummings et. al. [33], [34] applied the Volterra theory to X-31 airplane in order to identify transonic responses and pitching characteristics. Balajewicz et. al. [35] applied Multi-Input Volterra theory to transonic unsteady CFD simulations of a foil subject to heaving and pitching displacements K. Liu et. al. [36] applied Volterra Kernel identification to flapping 3D CFD simulations. However, only one maneuver for a small angle is identified and no information about the kernel is provided, so the application of the Volterra theory in the literature is very limited related to flapping wings.

Contribution
In fact, there are no analytical models that predict the effect of intermediate-range Reynolds and high induced angles of attack in ornithopters, in those ornithopters operating by large beating at low velocities [27], since hypotheses can be invalid. It is necessary to better analyze the vortex structures that are generated, in order to better understand the lift-propulsion generation mechanisms in the regime of large flapping ornithopters. Therefore, the identification of the Volterra model is proposed as a generalization of the classical linear models, without the need to separate between added mass and circulatory effects, and which can incorporate the effects that occur in the range of interest using empirical or numerical data. The main contribution of this study is the identification procedure of an unsteady aerodynamic ROM based on Volterra Theory for high-amplitude and low AR flapping wings at intermediate-range Reynolds. Moreover the analysis of the deficiency of the classical aerodynamics models at this conditions is also provided. Finally the validation methodology using Motion Capture data from flight experiments is proposed as an alternative to wind tunnel validation.

Contents
In section 2 the CFD scheme selected for the present study is presented, both model parameters, mesh and temporal discretization as well as the selected operating range for the aerodynamic database construction. Also a initial validation by literature data is provided. Section 3 provides an analysis of the unsteady 3D CFD aerodynamic database. Also presents an analysis of steady simulations in order to select the angle of attack range before unstable stall, and two-dimensional simulations at different heaving amplitude, in order to obtain the location of the surrogate airfoil. Finally, the underlying phenomena in the fluid is highlighted by vortex and pressure visualization. Then, the Volterra ROM selection, which depends on the kinematics of the surrogate airfoil, and the identification procedure is presented in section 4. The model is compared to classical Wagner, and the discrepancies are analyzed by frequency response. The identified kernel and verification of the identification is also provided. Finally, in the section 5 it is propose a validation methodology using Motion Capture System, and the validation for lift force is provided.

Computational Fluid Dynamics simulations
In order to gain insight into the phenomena of flapping flight, some authors proposed the resolution of Navier-Stokes equations. In reference [24] can be found most relevant CFD analysis of flapping airfoils, most of them focused on vortex visualization and correlation of two dimensional airfoils at low reduced frequency performing pitch and heave oscillations at different hinge points. However, the most challenging task on CFD simulation of 3D high rotation bodies is the efficient generation of an accurate mesh. The "remeshing" techniques were use to model the fluid volume elements along time [37]. After, Bos [37] among others proposed a mesh deformation technique based on radial basis functions to simulate the flow around an insect wing at Reynolds O(100) using a laminar model. Nowadays become relevant the use of local mesh adaptation or Overset techniques that joined to turbulence models provides an accurate tool to solve complex flows around flapping wings [38]. Unsteady 3D numerical simulations have been performed in the commercial software ANSYS Fluent R for the GRIFFIN prototype wing.

Simulations kinematics and operating range
It is widely known that the effects of chord based Strouhal and Reynolds aerodynamic numbers govern the fluid around flapping wings in forward flight, in addition to the geometric angle of attack. On the other hand, the amplitude of the flapping influences the aerodynamics if it has a considerable magnitude, that is why it is not considered in small perturbation's models, such as Theodorsen's model. At this study, the amplitude remains constant at a considerable magnitude. For that it will refer to the Strouhal as the chord-based Strouhal, defined below.
The wing is modelled as a rigid and thin ellipsoidal plate, which dimensions can be found at  The Overset Method has been used to simulate the displacement of the wing over the flapping, where two unstructured meshes are superimposed. The background mesh, provides the layout where slides the Overset mesh that contains the mesh of the wing surfaces. With this method, the meshes maintain their shape during simulation, avoiding the loss of quality produced by elements distortion that occurs in the dynamic mesh simulation. Figure 2a represents the dimension of the fluid domain for each of the meshes, which have been seen to be sufficient to cover the perturbations of the wing, taking into account the reduction of the distance to the walls with the movement of the wing. The wing is centered on the Overset mesh (2b), which is in the position shown in the figure at the initial instant. Note that the background mesh has been refined in the moving zone so that the displacement in a time step is less than the size of the cell for the critical case, that is, maximum frequency and speed. Mesh refinement has been increased without noticeable changes.  The simulation parameters are the fluid inlet velocity, the geometric angle of attack and the flapping frequency as Figure 3 presents. The mean angle of flapping or dihedral and the amplitude of the flapping of the wing has been fixed to 21.5 deg and 26.5 deg respectively. The first responds to stability needs, and others are chosen responding to bio-inspired design criteria. Simulations have been performed over the entire operating range of the ornithopter, which are shown in Table 2.  Reynolds-Average Navier-Stokes (RANS) model is chosen, precisely the k − ω with Shear Stress transport, a incompressible two equation Eddy-viscosity model commonly used in intermediate-range Reynolds aerodynamics. A pressure-velocity coupled scheme has been used, with a first order implicit time discretization for the transient simulation. Also one case has been simulated for the laminar model, obtaining negligible differences and an increase in computing time of more than double.

Database construction and verification
Once the operating range (α k , V k , f k ) and CFD model is selected the simulations are automatized based on the method proposed in Figure 4. The values of the simulation variables are provided by a Matlab R interface.  Figure 4: Aerodynamic database construction process. k is the simulation identifier (a combination of incident velocity, angle of attack, and flapping frequency), j is the time step identifier for each simulation.
For each simulation the timestep is selected automaticaly in order to obtain n points per flap, and N total flapping, hence f 1 = 1 f k n , where f k is the flapping frequency, n = 100 is the number of points per flap. The time vector will be j = 1 : N , where N = 3n. First, a steady state is simulated at each condition in order to initialize the flow before the unsteady simulation starts. For each combination of parameters, a good convergence of the temporal results after 3 flapping has been observed. This convergence has been improved by previous steady simulation of the fluid field. The fluid variables solution at each element is contained in R k . The steady simulation is iterated until the error became smaller than s = 1 · 10 −6 , obtaining a good convergence in the forces and moments at each time step. Then, the steady simulation is used as initial values to unsteady simulation. The kinematics of the wing is controlled by User Defined Functions (UDF), where a sinusoidal or chirp function (f 2 ) of the position of the wing Θ k j is enforced and the mesh is refreshed in each time step by this kinematics. The unsteady simulation iterates in each timestep until error is smaller than u . The total forces and moments are saved at each time step in inertial reference frame F k,I j and it is offline post processed in order to maintain a constant time step and disposed in aerodynamic reference frame (f 3 = gcd(dt k )), where gcd is the gratest common divisor. Also the effective angle of attack is calculated to complete de aerodynamic database. Equations f 4 , f 5 , f 6 are showed in detail in the next section. To finalize, the overall fluid variables R k j at each element are only saved at selected relevant points over the flap for vortex visualization in ensight-case-gold format due to the high computational cost.
The CFD scheme has been verified for a low-Reynolds high-amplitude critical case. The model used is an ellipsoid, which performs complex figure-eight movements. The results obtained shown in Figure 5 by the proposed Overset Method provide results in agreement with those offered by F. M. Bos in reference [37]. Black circles markers represents lift and drag coefficient from [37], black lines are coefficients from our CFD simulations reproducing the same conditions. The forces are in wing fixed reference frame.

CFD results: Flow visualization and surrogate airfoil
Approximately 70 simulations of flapping at constant frequency and 10 simulations varying the frequency in a linear chirp have been carried out in the previously defined range.
The effective or induced angle of a airfoil (see Figure 6) in the wing depends on the kinematics of the wing, and will be defined in the modelling section. The coefficients (Lift C L , Thrust C T and Pitch moment C M , aerodynamic coefficients C L = L

Infinite wing in heaving oscillations: Surrogate airfoil
Furthermore, two-dimensional simulations of the wing have been carried out, which will be used to find the equivalent airfoil position. The quality and characteristics of the mesh are equivalent, with a cell height of 2 · 10 −4 m approximately. While in a airfoil it has been shown that the aerodynamics depends mainly on the effective angle of attack at 3 4 of the chord [8], an analogous representative point can be found along the span. Hence the movement of a flapping three-dimensional wing has been reproduced by the simulations of the airfoils that compose it. This can be interpreted as the infinite wing in heaving equivalent to the three-dimensional finite wing. Figure 7 shows the amplitude of the lift coefficient of the simulated airfoils at different positions of the span ( , as a function of the Strouhal number, and the same results obtained for the three-dimensional flapping wing. the highest correlation is observed with the airfoil located at 70 % of the wingspan.

Finite wing in flapping
The aerodynamic coefficients obtained for a semi-wing, in "total aerodynamic" reference frame (in the direction of α ef f and its perpendicular ) for the last period of each simulation is shown, which has already stabilized. The data is presented versus the effective angle of attack, in function of the Strouhal (St) (from 0.1 to 0.84) and the static angle of attack (α) (from −8 to 20 [deg.]). Note that the coefficient of aerodynamic moment is measured at the center of the chord's root for convenience. Observing the lift results on Figure 8, an approximately linear section with the effective angle is observed, up to approximately 40 [deg.]. By increasing the Strouhal, the induced angle is increased, which has a limit of 90 [deg.]. The lift coefficient decreases/increases rapidly at at positive/negative large effective angles of attack values respectively. It is observed that the static angle of attack has little influence on the peak values, higher the Strouhal, since the angles induced by the wing are much greater than the static one. In all the simulations, a great temporal dependence of the effective angle has been observed, which justifies the use of memory-allowed models to identify them.
One flapping cicle Regarding drag coefficient (Figure 9), remains positive for most position of the wing, however, this does not mean that the wing does not produce propulsion, since the projection of these coefficients on traditional aerodynamic axes provides negative resistance in certain cases. The resistance is approximately zero for null induced angles, since it is a flat plate in which case the fluid would be parallel to it. It should be noted that if the forces are projected on traditional aerodynamic axes, the resistance would be different from zero for a null static angle of attack. In general, the drag depends on a quadratic form with the induced angle. The discrepancies that exist between the values for a positive and negative angle may be due to the symmetry break between the upstroke and downstroke offered by the dihedral of the wing. The pitch moment coefficient (positive head-up) is also non-linear ( Figure 10), mainly by an odd exponent. However, the temporal dependence is much less appreciable than in lift.

Drag coefficient [-] vs Effective AoA [-]
+ 0-1 ref. One flapping cicle, measured at root c/2 Regarding the flow, in Figure 11, it can be observed how at intermediate Strouhal, a Leading Edge Vortex detachment occurs, which influences the detachment of the Wingtip Vortex (WTV) creating complex structures downwards. While increasing Strouhal, also the trailing edge is detached, interfering with the same wing a flap after. It has been found that the LEV is detached at Strouhal higher than 0.3 at high angles of attack. Figure 12 shows the pressure contours during one stroke at different Strouhal number, 0.26 and 0.41 respectively. In the second case, the LEV detachment produces a low pressure gradient in the upper surface, providing an increase in lift.

Reduced Order Model identification
Since CFD is too computational expensive, a feasible approach widely used in aerodynamics is the development of a Reduced-Order Model (ROM). These techniques have been used for forces and moments prediction and control problems. A well-known time domain theory in the field is the Volterra Model. The control and design of ornithopters capable of performing maneuvers high precision models are needed in the operating condition range. The identification of reduced order models is proposed as a method to create models that fulfill those requirements, and which can be simulated to be validated in the proposed methodology.

Model selection
The Volterra theory of nonlinear systems is appropriate for aerodynamic ROM identification because it casual, time invariant fading memory characteristics. Its greatest potential lies in the ability to include system memory when computing outputs, as long as this memory is fading and weakly nonlinear [40]. Furthermore, it is well known in the field of fluid mechanics that previous states have an important role in aerodynamics in all its applications, revealed in the well known convolution integral, and Volterra series is a generalization of the convolution description of linear time-invariant operators to nonlinear operators [41]. While it has been used widely to identify aerodynamic models of fixed wing platforms from both experimental and numerical data, the application to flapping wing is actually limited [36], and it is necessary to gain insight into the application to aerodynamics of flapping wing identification regarding to states, kernels and parametrization. For a continuous time single input x(t) and output y(t) the model is presented as bi-truncated Volterra Series, which can be declared as a functional of the states history x: The limits of the integral in (1) indicate that the response of the system at instant t depends on the input variable in the whole previous time up to the instant t. This definition would be unapproachable for long simulations, and unrealistic because the effect of previous states on current loads fade. For example, the effect of a vortex detached on the aerodynamic characteristics of the wing in an instant t would diminish due to: First, the vortex fade away due to the viscous effect of the fluid, and second, this vortex is dragged downstream by the incident stream.
The Volterra series applied to aerodynamic data can be understood as a numerical generalization of classical aerodynamic theories, being able to take into account any phenomenon produced by current or past system states. The classical model of Wagner for a heaving input h and effective angle of attack α ef f =ḣ/V + α is shown below, where Φ is the Wagner's step response, where Φ(0) = 0.5 and Φ(∞) = 1.
As defined, the aerodynamic lift coefficient depends on all the past inputs, since the integral low limit is zero. However, the step response converges to stationary state after approximately τ = 20. Due to this, the condition can be relaxed to have a finite memory model using the Finite Impulse Response (FIR), shown in Equation (7) in discrete form,where M is the past time steps took into account, for the surrogate heaving wing, andᾱ ef f =α ef f b V is the non dimensional derivative of α ef f .

Identification
In the literature there are some techniques to identify the Volterra Models, which has been used in a variety of areas and applications. Regarding to identification method, the Cross-Correlation Method has been extensively used in signal processing and communications [42], [43], however the training set has different restrictions such as zero mean. The Impulse Unit Method become relevant in aerodynamics where Silva [31], [32] studied the application to unsteady aerodynamics of fixed wing platforms by experimental data. However the replication of an impulse change in AoA is challenging in CFD because of numerical instabilities. The identification process could be accelerate and optimized using basis functions in the field of aerodynamics, such as Orthonormal Basis Functions [44], however, a priori knowledge of the kernel it is needed. The direct identification has been used by Liu [36] to identify one maneuver of a flapping wing, where no information about kernels and method's details are provided, and it is chosen for the present method.
To test the ability of the proposed scheme, the Linear Volterra Model will be identified through the lift time series of the simulations. The derivative of the effective angle of attack in the position discussed above is chosen as the state, and the different coefficients of forces and moment as output. In the same kernel will be included all the phenomenon since the separation between circulatory and added mass effects has been shown arbitrary [30]. The model selected is shown in Equation 8, wherė H is the kernel to identify, analogous to the Finite Impulse Response (FIR).
By carrying out the necessary operations, the model become a linear system such as AH = B , which can be identified by various methods. If N is the length of the state and output (B) time series, the matrix A would have dimensions (N-M x M + 1), and it will multiply if the Nonlinear Volterra Model is chosen. However, nowadays computing power is rarely the main problem. In this case the main difficulty comes from the formulation of the identification method itself. This is badly conditioned by definition so the rank of the matrix is low. Therefore, direct identification using least squares is not feasible, and regularization techniques are needed.
To solve this, Ridge Regularization is proposed as a method to avoid over fitting due to multicollinearity. Although others have been used (Moore-Penrose Semi-inverse between them) this is the one that has led to the best results. The identification is defined as follows, where λ is the regularization parameter selected taking into account the lowest appreciable eigenvalue of A matrix, and I is the identity matrix of dimension (M+1).
After some tests, it has been obtained that the step response stabilizes after 600 time steps at dt = 1/1200, so the memory is set at approximately 0.5 seconds. The computation time for an identification set of a time serie with 24000 points is less than 5 seconds. The kernels have been identified for sets of simulations at different mean angles of attack and frequencies, keeping the upstream velocity constant. For each velocity, a kernel has been obtained, Figure 13.

Analysis
One of the main problems with the identification of the Volterra Model is the loss of physical meaning that occurs with the identification of the kernel. However, taking the frequency response of the Theodorsen model as a reference (Note that Wagner's model function Φ(τ ) in frequency domain corresponds to Theodorsen's model function C(s) ), the differences can be quantitatively analyzed and associated with aerodynamic effects. The Theodorsen's Model is shown in Equations 13 and 14 for α ef f andᾱ ef f as input respectively. Figures 14 and 15 show the frequency response of the data obtained from the CFD simulations, the Theodorsen model, and the Jones aspect ratio correction as a function of the effective angle and its derivative, wheres = ik is the non-dimensional Laplace variable.
C(s) = 0.5s 2 + 0.2808s + 0.01365 s 2 + 0.3455s + 0.01365 (15) The input for the 3D simulations has been chosen at 70 % of span, since it has been demonstrated that is the surrogate airfoil. Regarding the magnitude in Figure 14, it can be observed that for low and moderate reduced frequencies (k < 2), the classical models underestimate the lift force. This range is influenced, according to Theodorsen, by the circulatory forces in the airfoil, which are diminished by the effect of the aspect ratio. The WTV in flapping is generally higher than in heaving, so the force will be much less than predicted. In the high frequency zone (k > 3) an increase in lift is observed compared to that expected. This hyper-lift phenomenon has been pointed out by several authors as the cause of LEV detachment. Regarding to phase, in non viscid theory, a positive phase which tends to 90 deg is found in most reduced frequencies, showing that the influence of the first derivative is dominant, except for low frequencies. In the classical Theodorsen's model it can be observed that the forces due to the apparent mass gain relevance when the frequency of oscillation of the airfoil increases, producing a phase advance with respect to the effective angle of 90 degrees in the upper frequency limit, since it is directly proportional to airfoil acceleration. On the other hand, the magnitude is diminished for large frequencies, and the limit of the transfer function being 0.5 for large frequencies, in concordance to Wagner instantaneous lift. For small frequency values, the contribution of forces due to the circulatory characteristics of the airfoil becomes relevant.
However, a very remarkable fact is that the simulations tend to a lower value, around 30-45 degrees depending on the free stream velocity (The different lines seen in CFD correspond to simulations at different speeds, which when simulated up to a specific frequency , they reach different Strouhal.). This may mean that the influence of the angle and its first derivative remain equivalent. Taking into account the recent study by Taha [30], where it is shown how the circulation actually has a phase that tends to -45 degrees, as opposed to the controversially defined Theodorsen circulatory force that would tend to 0, it may be suspected that there are mechanisms not taken into account in the non-viscous theory that cause the phase to tend to a value less than 90 degrees under certain conditions. Due to this, it is not possible to identify a rational transfer function through our data as proposed by [25], since rational functions phase tends to a multiple of 90 deg, which multiplication factor depends on the order of numerator and denominator of the transfer function. Regarding the mean values of the aerodynamic forces that are not taken into account in the frequency analysis using the Bode diagram, they are shown as a function of the reduced frequency k, and the aerodynamic angle of attack α, from CFD Figures 16a and 16c, and from analytical theories Figures 16b and 16d. While peak values of the forces had a negligible dependence with the mean angle of attack, as much higher effective angles of attack were generated, in the means a great dependence with them is observed. The mean lift coefficient during a cycle is shown to be generally lower in CFD than that of the Wagner model, due to the WTV effect as explained before. Besides that, at small angles of attack and low frequencies the CFD result is negative. This may be because the vortexes generated during the cycle interfere with the circulatory effects that occur at small angles of attack. In general, both increase with the angle of attack, and decrease slightly with the Strouhal.
(c) Drag from Navier-Stokes simulation (3D). (d) Drag from Vortical Impulse Theory [12]. The average drag coefficient, or propulsion if negative, has been compared in the figure with that obtained by Feria [12] using the Vortical Impulse theory, for small induced angles and therefore small amplitudes. While this theory provides a monotonic relation with reduced frequency and little dependence on the mean angle of attack, the CFD simulations predicts a highly non-linear map, with dependence on the |α|. Also, the maximum thrust in the study range is at angles around 8 deg. and maximum reduced frequency, that is, no propulsion occurs for large average angles of attack. Also, in general, the propulsion predicted by CFD is lower.

Verification
The identified kernels have been verified using training chirp time series. The identified and CFD models are compared in frequency response in Figures 17 and 18, where Normalized Mean Squared Error is shown in Table 3.
It can be seen how the identified Volterra model correctly reproduces the frequency response of the system. The phase and magnitude is captured with great precision, however the amplitude at a very small Strouhal has a greater error. [  Table 3: NMSE of CFD and Volterra Model frequency response, in magnitude and angle. The maximum error is 6.6 % in magnitude and 3.4 % in phase. Note that Θ andΘ are the values from CFD and estimated from the identified model respectively.

Flight experiments: Validation.
In order to obtain estimates of the aerodynamic forces of the ornithopter, indoor open-loop flapping flights have been performed at different frequency, in the GRVC 2 facilities. Figure 19 represent a block diagram resume of the proposed approach to validate the identified Volterra Model by flight experiments.    The total forces and moments are measured by a precise Motion Capture System (MCS), which provides position and attitude with a frequency of 120 [Hz]. Redundant active markers fixed in the ornithopter are used, as shown in Figure 21 and the system is calibrated in such a way that a precision of 2 millimeter is reached at the position of the markers. The initial conditions for angle of attack, speed, and tail position both in pitch and roll, are such that the ornithopter is intended to start flight under trim conditions to promote sustained two-dimensional flight, and are calculated for each frequency manually. The flapping frequency remains constant at each experiment, and the velocity and angle of attack oscillates in the flight around the equilibrium position because of the absence of control. It has been proved that the flight is preformed sufficiently in a vertical plane after trim to assume that lateral forces can be neglected. To obtain the forces and moments, a two-dimensional movement is assumed, so the two-dimensional dynamic equations of flight are used. The inertia model of the ornithopter is shown in Table 4, wich has been calculated by high accurate CAD (Computer Aided Design) model. The measurement of the position of the markers is post-processed to obtain the attitude and force estimates through equations f 1−4 wich corresponds to Equations 16, 18, 17 and f 3 : V = √ u 2 + w 2 . The maneuver is simulated by the identified volterra model (f 8 ) using the measured effective angle of attack and velocity from flight (f 7 ) in order to compare aerodynamic forces.    The dynamic equations that allow to extract the lift, moment and drag (L, M, D) produced by the wing in position A previously described, are shown below. First the position measured in inertial reference frame (x h z h ) is derived and transformed to body reference frame (b). Then total forces and moment (F x , F z , M y ) in body frame are constructed by the body acceletarions and pitch (Equation 17), and tail influence is substracted in Equation 20.
For this, it is necessary to implement an aerodynamic model of the bioinspired tail (L t , D t ), which has been obtained through experiments in a tunnel of wind, and depends on the free stream velocity and relative angle of attack. More information about this model used in the scheme in (f 5 ) can be found in [2].
Once the data from flight is post-processed, the reconstructed model inputs, such as free stream speed, angle of attack and frequency are use to simulate the same maneuver by the aerodynamic model. The memory is adapted to flight time step which is ten times higher. In Figures 23a and 23 can be found various validation sets from flight compared to our model (Volterra Model) and the classical model of Theodorsen. The Volterra model reproduces well the amplitude and mean of the experiments, while Theodorsen overestimate lift amplitude and in lower level, mean lift. Validation set at k=0.9 Flight data Volterra Model Theodorsen (b) Reduced frequency 0.9.   Note that while the flapping frequency is selected a priori by flapping actuator commanded power, is measured in flight since depends on many variables like batteries level or environment. However, wing stroke position, which could provide the phase between the effective angle of attack and aerodynamic forces is not measured yet, so the phase is selected for best estimates.

Conclusions and future work
The aerodynamic forces and moments produced by an ornithopter elliptical wing in forward flight performing high amplitude flapping at intermediate-range Reynolds has been simulated by unsteady CFD computations. The results has been analyzed quantitatively and qualitatively using time series and graphical post-processing tools respectively. The Leading Edge Vortex is detached after Strouhal 0.3 is reached during middle down stroke.
Also it has been compared in frequency domain with the classical model of Theodorsen for a airfoil. For that, bi-dimensional CFD simulations have been performed where is shown that the three dimensional wing fits better with the airfoil located at approximately 70 percent of the semi span, performing heaving with the corresponding amplitude. So regarding to modelling, the effective angle of attack located at this position span wise and three quarters chord in chord-wise direction is chosen as input for the aerodynamic model. It has been found that the magnitude is overestimated by classical models at low reduced frequencies. The phase tends to approximately 30 degrees instead of 90 degrees of the Theodorsen's model.
The CFD simulations has been identified by a Volterra Model using the linear kernel, where a memory of 0.5 second has been found sufficient. A ridge regression has been found the best method.
Finally, the identified model is validated through aerodynamic forces reconstructed from flight data positions and attitude measures by a high fidelity Motion Capture System.
The model is capable to be used to simulate the aerodynamic forces acting on an ellipsoidal wing ornithopter performing high amplitude flapping at high angles of attack and intermediate-range Reynolds number, for any time serie of effective angle of attack (Including flapping frequency, incident velocity and mean angle of attack.). The computation time of the model is negligible, so it could be used to online estimation. In future, it would be measured the position of the wing in flight to validate the phase between effective angle of attack and the aerodynamic forces and moments. This could be done using MCS and placing markers on the wing to become an skeleton instead of rigid body. Also, a lateral forces Volterra Model could be identified, a nonlinear extension of the Volterra model could be done, or the append of elastic terms to the aerodynamic model. Regarding to future work, a similar conclusion than other ornithopter researchers has been reached, which is that the state of the art of ornithopters aerodynamics and flight mechanisms is lower than the insects robots. So there are various topics which still needs to be studied, like the effect of the lateral aerodynamic in flight dynamics, optimum kinematics of birds wings among other.