Optimal Elastic Wing for Flapping-Wing Robots Through Passive Morphing

Flapping wing robots show promise as platforms for safe and efficient flight in near-human operations, thanks to their ability to agile maneuver or perch at a low Reynolds number. The growing trend in the automatization of these robots has to go hand in hand with an increase in the payload capacity. This work provides a new passive morphing wing prototype to increase the payload of this type of UAV. The prototype is based on a biased elastic joint and the holistic research also includes the modelling, simulation and optimization scheme, thus allowing to adapt the prototype for any flapping wing robot. This model has been validated through flight experiments on the available platform, and it has also been demonstrated that the morphing prototype can increase the lift of the robot under study by up to 16% in real flight and 10% of estimated consumption reduction.


I. INTRODUCTION
F LAPPING wing UAVs are a very active line of research due to its great potential for application in natural or nearhuman environments, thanks to their safety and low environmental impact (visual and noise). The bioinspired operation of these platforms, as it is the case of perching, requires a high level of automation, sensors and soft devices onboard, for which a large payload capacity is needed [1], [2]. On the other hand, birds modify the shape of their wings to maximize lift either in low-speed flight or when carrying an extra weight, inspiring researchers to design and analyze new morphing systems. Morphing systems can be classified into active and passive. While the first invest an extra amount of energy in the deformation of the wing, passive systems are deformed via the fluid-solid interaction and natural movement of the UAV. In the case of flapping-wing robots, the design of such devices are less explored due to the complexity of their aerodynamics, such as high unsteady and three dimensional viscous flows, in addition to the large inertial forces that they must withstand during the flapping at usual frequencies. The current state of the art on ornithopters Manuscript  wing morphing is extremely narrow, including mainly active morphing such as rigid-bar mechanisms as in [3] and [4] or more complex that mimic bats kinematics [5] or use active morphing by shape memory alloys (SMA) [6]. Regarding passive morphing, the tendency is to introduce an asymmetry to compliant mechanisms in order to decrease the projected surface of the wing during upstroke, and therefore the negative lift, leading to an increase in the wing mean lift force. However, the proposals are very limited in the ornithopter literature, and void when flight validation is considered. Although there are science study cases for sweeping in MAV [3], most research, and this one, identify bending as the most evident deformation in the natural flight of birds. Billingsley [7] first test the implementation of a torsional spring in an ornithopter wing for bending. A stop mechanism maintained undeformed the wing during downstroke for lift maximization. The spring elasticity was selected experimentally and qualitatively. The conclusion of the study in bench, is that the mechanism makes the average lift increase noticeably, however the propulsion decreases in such way that cannot generate forward velocity and hence, a sustained flight was not possible. In the conclusions, the authors highlighted the need of a 2 DoF dynamics model the of the wing for sizing and design of the spring, which is one of the contributions of the present work. At the same time Mueller [10] applied it for the decrease of the flight velocity and therefore the increase of the maneuverability of Micro Aerial Vehicle (MAV), with a validation on a bench setup. Subsequently, in the line of smart material and structures, Wissa [8] developed a mechanically complex compliant spine, in order to improve the steady level flight performance. The spine was designed in order to bend in upstroke while minimizing the mechanical stress, and it was placed at 39% of the span. An increase in mean lift was observed in the bench tests, however, no extrapolation to real flight has been found. Later in [9] the authors test the spined wing prototype in tethered (constrained) flight, which data obtained through a VICON system is used to analyze the power requirements and acceleration of the center of gravity. However, no lift data were reported. Moreover, in future work the authors highlighted the need of a quantitative comparison between the free and constrained flight test results.
In summary, the literature shows how few preliminary tests of wings with passive morphing have been carried out, which, however, have not led to stable free flights. The proposed passive morphing system compared to the actual state of the art is summarized in Table I. Notice the increase of lift by 16% validated with flight data and the decrease of consumption of 10% estimated with the model developed and validated.
In this work, the implementation of a free elastic joint is proposed, in which a bias angle performs an asymmetry between downstroke and upstroke, thus (passively) modifying the projected area and hence the generated lift smoothly. In addition, to avoid manufacturing errors highlighted in all the references above, a rigid link has been introduced within the compliant itself to restrict movement only to the degree of freedom of interest and stabilize the wing beat. An aeroelastic model of the wing has been developed and experimentally validated. The model yields the expected aerodynamic loads for a combination of the design parameters. Moreover, it allows us to automate the numerical optimization of the mechanism design and sizing to obtain the best lift performance, thus avoiding any manual or heuristic procedure to select the 'optimal' design. This optimization is posed as maximization of the objective functional given by the total mean lift over a flapping period, subject to flight regime and dynamic constraints. The selected parameters are the elasticity and position along span of the joint and the bias angle.
The contributions are enumerated in detail below: C1. A novel aeroelastic model of a flapping wing with an elastic elbow and its experimental validation throughout flight data. The parameters are intrinsic to the elastic joint model and then can be extrapolated to any ornithopter even at low Reynolds and large flapping amplitudes. C2. Design through a constrained optimization of the passive morphing wing, accounting the spring elasticity, location along the wing and bias spring angle. C3. The mechanism is validated in flight. Quantitative improvement of mean lift is reported based on the experimental flight data in a large scale flapping wing UAV.
To the best of the authors' knowledge, no contributions of type C1 and C3 for passive morphing in large scale flapping wing UAV have yet been reported. The E-Flap bird-scale flapping wing robot [1] (see also related work [11], [12], [13]), is the platform used for testing the proposed morphing mechanism, presented in Fig. 1. With a total weight of 0.650 kg, large payload (up to 100%), large amplitude flapping wing 60 deg and a 1.5 m of wingspan.
The paper is structured as follows. In Section II, the novel aeroelastic model is derived for a general flapping elastic wing with compliant biased joints, and preliminary results are analyzed. In Section III an optimization algorithm is proposed with the obtained model. The optimal parameters for the E-flap platform are used as design layout for the experimental validation using a Motion Capture System, which is presented in Section IV. Finally, conclusions are included in Section V.

II. MODELLING AND SIMULATION
Performance of elastic flapping wings is a purely dynamic aeroelastic phenomenon. The need for modeling and optimization of this system arises from the inherent difficulty due to the non-linearities produced by both non-stationary aerodynamics and the complexity of the system, making its optimization by trial and error either inaccurate or unapproachable. Nowadays there are numerical aeroelastic simulation software available by combining CFD and FEM in FSI schemes, however the high computational load that they demand makes them practically infeasible for optimization tasks in unsteady flows. On the other hand, existing aerodynamic and elastic models are of medium reliability when used for numerical optimization without excessive computational demand. Thus, in this work and to solve all the aforementioned drawbacks, a novel model is proposed that combines a simplified elastic model with a reduced-order aerodynamic model validated in flight. The aerodynamic model selected is the Volterra Model developed in previous research in [13]. As it was thoroughly explained in there, Volterra's model is causal, time-invariant and includes fading memory characteristics. Thus, it is able to reproduce the unsteady aerodynamics [14] produced by the large flapping amplitudes of the ornithopter, with accurate CFD simulations while maintaining the physical sense of the identified kernel.
The elastic flapping wing is symmetric in yz plane, so only the half wing is modeled by two rigid bodies, inner and outer as shown in Fig. 2, whose total uniformly distributed mass is M = mD where m is the mass density in kg/m and D the half span. The inertia of the inner body (wing section 1 in Fig. 2) respect to flapping rotation edge becomes I 0 1 = 1 3 md 3 , and for the outer body (wing section 2 in Fig. 2 3 with respect to its gravity center, where d denotes the span of the inner part and hence the location of the joint. Both sections are assumed to be square plan shape for simplicity, however, any wing inertia can be implemented without loss of generality, In fact, a thorough analysis performed resulted that changes in the planform wing slightly affect the optimal elasticity. The inner body is considered joined to an inertial frame by a frictionless 1 DoF joint, and linked to the outer part by a torsion spring of elastic constant K, as shown in Fig. 2. The air conditions are defined by the inlet velocity V and the geometrical angle of attack α 0 . The variables considered for the model are (see Fig. 2): the flapping angle θ 1 (t), which is the angle of the inner body with respect to the inertial frame; the elastic deflection wing degree of freedom (DoF) ψ(t), selected by convenience as the angle of the outer body with respect to inertial frame. The latter can be decomposed into three contributions: ψ(t) = θ 1 (t) + θ 2 (t) + φ, where θ 2 (t) is the elastic deflection and φ the bias spring angle or the zero-load deflection, which is the cause of the asymmetry between upstroke and downstroke.
The system is exposed to the unsteady aerodynamic and inertial forces, F A,i and F G,i , i = {1, 2}, and the elastic and motor moments, M S and M F , as shown in Fig. 2. The aerodynamic force is applied at 3/4 span of the body [13], while the gravity force is applied at 1/2 span. The drag force is not included in the aerodynamic force as a trade-off between the model complexity and accuracy achieved. Indeed, by assuming a common efficiency L i /D i = 10 and worst case scenario α = 30 deg, the vertical force taking into account the drag becomes that estimates a maximum error of 5%. The lift force, is modeled for each wing section separately by an unsteady aerodynamic model [13], which upon adapting the equations the lift leads to where q i is the dynamic pressure, q 1 = 1/2ρV 2 cd and q 2 = 1/2ρV 2 c(D − d), ρ is the air density, c is the mean wing's chord and Φ(t) is kernel function of Volterra model, which depends on the non-dimensional time defined as t = t 2 V /D and on the induced angle of attack. As common in unsteady aerodynamics [14], the left part of (1) refers to the instantaneous lift and the right part is the 'memory' term in integral form, that is the lift contribution of the past states. The induced AoA, α 1 and α 2 , are composed by the classical geometric angle of attack α 0 and the one induced by the vertical component of the flapping displacement at 3/4 of the span in each section. For clarification regarding induced angle of attack see [13]. By using the relative motion kinematics, and defining the angular velocities in each section asθ 1 andψ, the induced angles of attack are given by The dynamical model is derived in the well-known Lagrange framework with generalized coordinates q = (θ 1 , ψ) as with the Lagrangian defined as the difference between kinetic and potential energies and Q being the generalized load moment associated to each DoF. Note that even though θ 1 is imposed, the equation associated to it is used to calculate the required flapping moment M F , and subsequently the power. Let us consider the following non-dimensional parameters and momenta: mass μ = 4 m/πρD 2 , aspect ratio AR = D/c, gravityḡ = gD/2 V 2 , joint positiond = d/D with ξ = 1 −d, and the elasticity and flapping momenta as Γ = K/ 1 8 ρV 2 πD 3 andM F = M F / 1 8 ρV 2 πD 3 , respectively. Accordingly, the nondimensional coordinates areψ =ψD/2 V andψ =ψD 2 /4 V 2 , and analogously forθ 1 . Thus, taking the discrete time ast = jdt for j = {1, . . ., N}, with dt the time step, the resultant nondimensional dynamical model becomes where the following functions have been defined the memory-like functions as and the constants C i , i = {1, . . ., 10}, are defined below Finally, the flapping kinematic is enforced as θ 1 (t) = A m + A f sin (2πf t), where A m and A f are the dihedral and flapping amplitude, and f is the flapping frequency. The system has been numerically simulated by finite differences according to the Algorithm 1, where N is the number of timesteps and p 1 = 2, p 2 = 3 are the orders. Note that the term related to k = j in (6) has been extracted, so that the memory-like term (6) can be calculated prior solving the equation (see Algorithm 1). A sensitivity analysis have been performed, obtaining an optimal timestep of dt =0.01 s. Besides, the aerodynamic mean lift is stabilized after 3 cicles aproximately, so the 4th cycle is used to calculate the mean lift. Importantly, with this computational approach, the simulation of 1 s of real-time flapping is performed in 0.25 s of simulation time on a regular laptop, making then onboard simulation feasible.

III. OPTIMIZATION AND ANALYSIS
In this section, an optimal design of the passive morphing wing for the actual platform E-flap [1] is developed. To this end, the optimal solution is obtained through the aeroelastic model of the wing with joint (3)-(4), which is validated with flight data. As will be shown, a suboptimal design causes the joint to decrease in performance. Moreover, by means of the Algorithm 1: Aeroelastic Model Simulation.
for j = 1 : N dȱ From eq. (1) end for flight data obtained in [13], a thorough comparison with a wing with no joint of [13] is also carried out. The underlying idea for the optimization is to design wings with elastic joints for the E-flap prototype to improve its payload. Therefore, there are fixed parameters as: half wing mass m = 0.1 kg, mean chord c = 0.36 m and half span D = 0.75 m. Table II collects all the 'tunable' parameters for the optimization problem. However, to reduce the computational burden, it is reasonable to fix also the regime with which the E-flap was designed. Hence the last 5 of Table II have been set to: V = 4 m/s, α 0 = 20 deg and gait A m = 21.5 deg, A f = 26.5 deg and f = 3 Hz. Notice that, the gait (A f , A m ) obviously will have influence, but our preliminary analysis shows that is not major. Thus, the design of the wing with joint has been optimized using the three parameters with the greatest influence, which are: joint elasticity K, bias spring angle φ and joint position on the span d. To that end, let Θ = {K, φ, d} be the outcome vector and denote the objective function to be optimize the total mean lift as L m = 2(L m 1 + L m 2 ), L m i > 0, from (1). Note that we denote the mean value over time with superscript m, so L m 1 and L m 2 correspond to the mean values of (1). The constrained optimization problem is stated as follows max Θ {L m (Θ,q, q))} subject to Remark 1: It is worth noting that, at first glance, it seems that the objective L m is an affine function in the dimensionlessd, but it appears in some coefficients C i . In fact, it is not possible to establish analytically any convex/concave property, because the need of the integral curves of (3) and (4). Additionally, this fact precludes the use of numerical approaches based on the vertexes of Θ. However, as it can be foreseen by physical considerations, the numerical solution obtained clearly shows an optimum that is experimentally corroborated.
Optimal design analysis. In Fig. 3 the mean lift (L m ) and required power (P m = mean(M Fθ1 )) are shown, as function of the main design parameters Θ. For the sake of comparison, the results obtained with a rigid wing, i.e. no elastic joint, have been also included. These results differ from the rigid wing results in [13] which were obtained for φ = 0 only, while here φ varies throughout the search range. As expected, for a constant joint position d, there is an optimal combination of elasticity K and bias spring angle φ that provides higher mean lift than the rigid wing design. This optimum (shown as a white dot in Fig. 3) tends to diminish elasticity as joint position increases, in order to compensate the decrease of inertia on the wing section 2 (see Fig. 2), maintaining a similar gait. Regarding required mean power, it can be found a minimum d(%) ≈ 15% from which the required power is smaller than the one of the rigid wing, that is, providing a more efficient elastic wing. The optimal design parameters at the design point are the elastic constant K = 3.49 Nm/rad, bias spring angle φ = −1.38 rad and joint distance d = 0.34 m. With this design, the model predict and increase of lift up to 33% (represented by circle marks) with respect to the rigid wing, 1.16 kg and 0.87 kg, respectively. Furthermore, the power consumption is reduced by 20% with respect the rigid wing, 39.18 W and 49.22 W, respectively. However, due to instabilities reported in experiments, the bias angle has been limited to 60 deg. The optimization has been repeated with this restriction yielding the actual design at K = 6.18 Nm/rad, φ = 1.04 rad and d = 0.3 m, thus providing an increase of 28% in lift and 10% decrease in consumption estimated. Note that there are other more efficient designs in terms of power, however, in the current design mean lift has prevailed over consumption, motivated by the increase of the payload.
Remark 2: An interesting observation is that the optimal configuration exhibits a gait similar to the one of birds in their flights. To show this, the lift and demanded torque over the flapping and the deformation of the wing φ are shown in the Fig. 4 , for the optimal configuration. Note that it is 'optimal' to  keep the wing horizontal, i.e. deflection angle φ ≈ 0, during the downstroke while maximizing the angle on the upstroke.
Analysis of performance. In order to prevent re-design due to abrupt changes in the elastic wing performance when operating conditions vary during flight, the performance of the actual design has been analyzed through the flight envelope (V, α 0 , f). Thus, the wing performance over the envelope is shown in Fig. 5. The lift increases primarily with angle of attack and velocity, having little influence on the frequency. In contrast, required power, which is related to consumption, is highly affected by frequency. The increase of the frequency by 1 Hz (33%) duplicates the required power, so it should be taken into account in the flight control system.

IV. EXPERIMENTAL VALIDATION WITH FLIGHT DATA
Prototype. To demonstrate the validity of the optimal design and model proposed, an elastic wing prototype has been manufactured for the E-flap platform with the aforementioned the optimal parameters. The CAD design for the joint is shown in Fig. 6 and it is composed by two pieces manufactured by 3D printer PLA, designed to hold in the desired angle a double torsional spring made of alloy steel. Furthermore, a customized aluminum rod has been inserted into the main wing stringer, forcing the deformation to occur in a single axis and minimizing the instabilities reported in the literature. The weight of the joints is approximately 0.05 kg, resulting in a total weight of 0.7 kg of the UAV, including all avionics.  Experimental setup. In order to validate the proposed model, a set of flight experiments has been carried out. A motion capture system has been used, with 28 infrared cameras distributed in an indoor space of 20 × 15 × 7 meters. The system has been calibrated in such a way that the position of each marker is obtained with an accuracy of less than 1 mm. A total of 14 passive markers have been located through the ornithopter, measuring the position and attitude of 5 bodies, the fuselage and the two sections of each half wing, as shown in Fig. 7. The strategy followed when positioning the markers has been to minimize the number while maintaining the precision of the desired measurements. The wing sections 2 are represented with only 3 collinear, however it does not affect the accuracy since the only relevant position is the rotations of the bodies with respect to their axes x r 2 and x l 2 . Thus, the markers position represent an unambiguous body. The autopilot used to maintain longitudinal stable flight is the one presented in [1] and [13] (see also [11]). The attitude and position measured in flight are used in the autopilot, via local network, and also to reconstruct the aerodynamic forces offline. On the other hand, the measurements of the position of the wing bodies are used to calculate offline the flapping angle θ 1 (t) and ψ(t) described in section II, to validate the proposed model. Attached to the letter a video is included, to show the qualitative difference between the flight with and without optimized wing (see also at https://youtu.be/kwuW8cfy-MI).
Experimental results. More than 30 flights have been performed, of which the 7 best flights for validation have been selected, at different frequencies. In Fig. 8 snapshots over a   Table III.   Table III.  Fig. 9, where x h , y h , z h is the inertial reference frame. The main body variables of a typical flight are shown in Fig. 10. In top figure, the linear velocities in body frame u, v, w and the velocity magnitude V = √ u 2 + v 2 + w 2 . As can be observed, after the launch, the UAV is stabilized in a operating condition. These variables has been calculated by the derivation of the center of gravity marker position of the body. In bottom figure, the body attitude in Euler reference frame, roll, pitch and yaw (Ω x , Ω y , Ω z ) and the geometric angle of attack α 0 = arctan(w/u). The total lift has been reconstructed by the 2D flight dynamics from the total forces, by removing the tail contribution assuming a tail model thoroughly described in [12]. The validation methodology is synthesized Fig. 11. Once reconstructed the timeseries from flight (f, V, α 0 ), then the lift and wing deformation together with the model described in Section II is simulated for the flight conditions and the actual joint design configuration. Finally, the simulated lift and deformation is compared to the experimental data, as shown in scheme on Fig. 11. For visualization, the result of one flight is shown in Fig. 12. The deformation of the wing, the wing gait, and the mean lift are simulated with accuracy. The simulation overestimates the lift amplitude, however this is not a critical fact in the payload estimation, which mainly depends on mean lift. The results are presented in Table III. In the left side, the flight characterization of the experiments with the elastic wing, flapping frequency f and velocity V , are presented beside the reconstructed mean lift force L m exp . The operation mean frequency increases from 3.08 to 3.83 Hz, which increases mean flight velocity. Then, the same flight has been simulated with the aeroelastic model described in Section II for validation, and the mean lift L m sim can be found in the table. The accuracy of the aeroelastic model can be predicted by the RMSE between the simulated and experimental lifts, obtaining an average error value of 11.9%. Besides, in Table III, the flight data of the prototype with the rigid wing is also presented. The flights have been selected from [13] so that the frequency is similar to that used with the elastic wing, to ease their comparison. As it can be seen, the performance of the rigid wing makes the velocities higher. The last column presents the relative comparison of lifts, such that the elastic wing provides 16% higher lift than the rigid one, on average. This is similar to the one reported by Wissa [9] in bench, however the present mechanism is simpler and the optimization takes some seconds in contrast to the computational demand FEM. In addition, rigid wing flights have a higher average velocity, which means that the lift coefficient is considerably smaller than the one of the elastic  III  ELASTIC WING EXPERIMENTAL AND SIMULATED PERFORMANCE COMPARED TO RIGID WING EXPERIMENTAL PERFORMANCE wing, as presented in Fig. 13. Taking this into account, it can be deduced that the elastic wing at the same velocity will produce an even greater lift increase, closer to the simulation predicted 28%. This may be due to the fact that the elastic wing decreases the thrust with respect to the rigid one, reaching lower velocities in flight, which could be resolved with a decrease in aerodynamic drag.

V. CONCLUSIONS AND FUTURE WORK
In this letter, we provide a functionality upgrade of flapping wing UAVs to increase their payload capacity, including a novel bioinspired flapping-wing high-lift device, its aeroelastic model for design optimization and its validation with experimental flight data. In addition, the optimized elbow distance result is found as an avian trend in nature. The proposed aeroelastic model serves to predict the wing deformation, lift and power required throughout the flapping. This model is based on a ROM from CFD and an elastic model for the joint, that can be used to obtain an optimal design for any ornithopter, even operating at high angles of attack and low Reynolds. The optimal design yields a smooth gait and the instabilities reported in the literature have been minimized by design. The wing has been characterized by indoor flight experiments. Measurements of onboad markers by a precise motion capture system are used to reconstruct deflection and lift, validating the proposed model with a 10% error in the mean lift. Using flight data from rigid wing experiments and novel flexible wing experiments, the total increase in lift is found to be 16%. To the best of our knowledge, this is the first in-flight (no tethered) validation of an elastic elbow in a bird-scale flapping wing robot. Future work includes an onboard power meter to estimate the real efficiency increase in flight.