Nikola N Gavrilović

Emmanuel Bénard

Philippe Pastor

Jean-Marc Moschetta

Performance Improvement of Small UAVs Through Energy-Harvesting Within Atmospheric Gusts

Fixed-wing mini aerial vehicles usually fly at low altitudes often exposed to turbulent environments. Gust soaring is a flight technique of energy harvesting in such a complex and stochastic domain. Presented work shows the feasibility and benefits of exploiting non-stationary environment for a small UAV. Longitudinal dynamics trajectory is derived showing significant benefits in extended flight with sinusoidal wind profile. Optimization strategy for active control has been performed with the aim of obtaining most effective set of gains for energy retrieval. Moreover, three-dimensional multi-point model confirmed feasibility of energy harvesting in a more complex spatial wind field. Influence of unsteady aerodynamics is determined on overall energy gain along the flight path with active proportional control. Most contributing aerodynamic parameters are identified and suggested as basic objective function of an UAV design for energy harvesting in gusty environment. In addition, passive approach of control related to structural dynamics is investigated, pointing out its potential and possible improvements with aero-elastic tailoring.

I. Introduction

R ecent works driven by the unfailing experience from nature has shown significant amount of energy available in atmosphere. This energy is coming through the form of vertical air motions described as spatial gradients such as thermals, shear layers, orographic lift, and short-period temporal gradients as for example gusts. There are strong indications that birds use their feathers for sensing of flow perturbations over their wingspan. [START_REF] Scott | Bird Feathers -A Guide to North American Species[END_REF] Being fluffy and subjected to uttering provoked by small disturbances, birds have natural sensory system which enables them to feel flow disorders along wing, even before inertial response. Eventual immediate action due to surface pressure fluctuations by modifying wing geometry or profile curvature allows quick and effective response in gusty environment. However, for a variety of reasons, it is understood that identical copies from nature to man-made technologies are not feasible. Instead, a creative inspiration and conversion into technology is often based on various steps of abstraction.

One of the key objectives in the process of UAV design is long endurance flight. The flight technique inspired from albatrosses in form of optimal trajectory in presence of wind gradients promise neutral energy cycles. On the other hand, gust soaring promise quiet attractive flight strategy for endurance enhancement as it does not require any path planning or veer of the course. The performance of small UAVs being constrained by on-board energy due to their limited size can be significantly enhanced by specific flight strategies according to expected atmospheric formations or even continuous disturbances. UAVs have already proven the useful side of their exploitation and our goal here is to extended their operability without disturbing defined mission objectives in any flight case scenario. Most of the energy harvesting methods rely on active control system that detect and exploit the energy of atmospheric turbulence through intentional manoeuvring of the aircraft. It was demonstrated that exploitation of the energy available in atmosphere requires responses that are usually beyond human capabilities, implying autopilot with active control. Despite the fact that conventional approach was based on mitigation of gust effects, such exploitation presents a considerable opportunity to improve all ready degraded performance of small UAVs.

This paper aims at providing a clear picture of the physics involved in the energy extraction mechanisms in the case of gust soaring, implementing three-dimensional and unsteady effects. It takes advantage of nowadays powerful computational resources approaching to more realistic flight case scenario and evaluation of possible benefits. The following sections of the paper are organized as follows. The section II covers recent developments in the filed of soaring strategies related to gusts. Section III briefly describes artificial generation of gust environment pointing out benefits of empirical methods. Section IV exposes system of differential equations of motion, active control strategy and dynamic energy transfer between aircraft and atmosphere for unsteady aerodynamics. Section V describes the methodology of modelling the structural dynamics, pointing out most important structural parameters that contribute to energy gain. Section VI quantifies the benefits of gust soaring detecting aerodynamical parameters of the aircraft which affect energy gain. Finally, conclusions are presented in section VII.

II. Related Research

P revious work performed by Patel and Kroo 3 examines the way of energy harvesting from vertical tur- bulence by active control of flaperons. Work from Langelaan 4 and Depenbusch 5 uses longitudinal vehicle model with optimized control in both sinusoidal and Dryden turbulence profile. The development of energy equation by Lawrence 2 and Langelaan [START_REF] Langelaan | Biologically Inspired Flight Techniques for Small and Micro Unmanned Aerial Vehicles[END_REF] showed that certain manoeuvres are necessary in presence of wind derivatives in order to extract energy. Optimal trajectories for energy extraction in form of neutral energy cycles in presence of wind gradients are investigated by Zhao [START_REF] Zhao | Optimal Patterns of Glider Dynamic Soaring[END_REF] and more recently Bonnin. [START_REF] Bonnin | Energy-Harvesting Mechanisms for UAV Flight by Dynamic Soaring[END_REF] Exceptionally complicated urban spaces pose a challenge for autonomous operation of UAV. In the urban environment the characteristic scale is on the order of few meters, which is due to a complex interactions between buildings, trees and living world. Simulation tool that captures unsteady aerodynamics of flight through such a complex terrain has been presented and demonstrated by Cybyk [START_REF] Cybyk | Unsteady Airflows and Their Impact on Small Unmanned Air Systems in Urban Environments[END_REF] . All strategies previously described omit advance knowledge or prediction of wind velocity field ahead of UAV. However, method for sensing flow disturbances in front of mini UAVs and using the output signal for further control has been demonstrated by Mohamed. [START_REF] Abdulrahim | Developement and Flight Testing of a Turbulence Mitigation System for Micro Air Vehicles[END_REF][START_REF] Watkins | A Feasibility Study of Bio-Inspired Wing-Surface Pressure Sensing for Attitude Control of Micro Aerial Vehicles[END_REF] The control strategy has been developed for roll axis as most sensitive to wind turbulence. Recent experiments by Watkins 10 considering measurements of wind components on locations along span-wise confirmed statement about spatial variation of turbulence magnitude which is in correlation with actual needs. Being three-dimensional, turbulence scales larger than wingspan would result in only pitching attitude of the aircraft. However, case of turbulence smaller than wingspan leads to unequal lift distribution and need for control of oncoming roll and yaw moment. Method for atmospheric disturbance effect on roll and yaw motions has been proposed by Ringnes 11 based on well-known Prandtl lifting line theory. Moreover, beside already described methods of gust energy extraction using active control, some benefits could be achieved using passive approach that makes the use of longitudinal stability of the aircraft and the dynamic response of the structure, particularly of the wing by Ironside [START_REF] Ironside | Modelling of Wing Drag Reduction due to Structural Dynamics in Atmospheric Gusts[END_REF][START_REF] Ironside | Simplified Modelling of wing Drag Reduction due to Structural Dynamics and Atmospheric Gusts[END_REF] and Mai. [START_REF] Ulv | The effect of aeroelasticity upon energy retieval of a sailplane penetrating a gust[END_REF]

III. Wind model

A ll the disturbances present on Earth or even far away can contribute to wind formations. Either unequal heating of the ground provoked by the clouds or natural obstacles as cliffs, mountains and valleys or even sea currents are provoking uneven pressure field that drive the motion of air. Beside low altitude, UAVs usually fly in urban environment surrounded by buildings, trees and other obstacles, where the characteristic scale is on the order of few meters. Such a complex surrounding implies intricate interaction between terrain geometry, physical conditions and varying meteorology. It is well known that turbulence intensity increases nearing the ground strongly influenced by the terrain, thus changing the conditions comparable with those at high altitudes. From field experiments it is well known that undisturbed wind velocity is variable in space, time and direction. A model of the turbulent wind field suitable for calculations requires good representation of both the temporal and spatial structure of turbulence. The most adequate method to simulate a turublent wind field would be to solve Navier-Stokes equations of an atmospheric flow bounded from below by an aerodynamically rough surface. However, the computational cost would be enormous. Alternative could be Large Eddy Simulations (LES) as an approximate solution to the Navie-Stokes equations, where the smallest scales are not solved directly but modelled. Still even simplified alternative requires big computational power. Therefore, empirical description is generally used for turbulent flight using spectral and coherence description of turbulence.

General algorithm to simulate three-dimensional field of three components of wind velocity fluctuations is developed by Mann. [START_REF] Mann | Wind field simulation[END_REF][START_REF] Mann | The spectral velocity tensor in moderately complex terrain[END_REF] The method builds a model of spectral tensor for representation of atmospheric turbulent boundary layer. It leads to natural representation of three-dimensional turbulent flow of hight cost effectiveness compared to alternative Large Eddy simulation or Navier-Stokes equations. The spectrum used in simulation is presented by Kaimal 17 where α = w x , w y , w z and A, B are adjustable constants that depend on the chosen length scales. Method for generation of a single wind time series from a Kaimal spectrum is proposed by Branlard. [START_REF] Branlard | Generation of time series from a spectrum: Generation of Wind time series from the Kaimal spectrum, Generation of Wave time series from Jonswap spectrum[END_REF] The generated wind field provides three-dimensional grid of desired spatial resolution which will be used later as a environment of flight for gust energy extraction. The characterisitcs of generated profile are compared with available database on wind characteristics 20 which provides overall information of wind measured at various mast heights near the wind turbine stations at many places in the world. Reference locations for turbulence energy were Tjare spot in Denmark described as pastoral, flat landscape and San Gorgonio in US with rolling hills covered with bushes and small trees. The corresponding referent wind profiles were taken at the mast height of 30 and 40 meters of altitude. These profiles are the reference for parameters of Kaimal spectral formulation in order to generate wind field as closer as possible to reality.

F α (k 1 ) = Aσ 2 u 1/k 1m 1 + B(k 1 /k 1m ) 5/3 (1)
The key parameters associated with external flows that drive urban environment include the following parameters according to Cybyk: 12 intensity of turbulence, turbulence length scales, surface roughness, Reynolds number and Richardson number.

IV. Flight dynamics

H ere we expose the system of differential equations that describes the physics of flight in gusty environ- ment. The first derivation of equations is done for the longitudinal plane. Also, some efforts are exposed in order to reveal dynamic energy transfer between the aircraft and wind. Firstly, the simulations are performed with sinusoidal vertical velocity wind profile for various frequencies taking into account unsteady behaviour of the aircraft. Secondly, more realistic multi-point model has been constructed for flight within three-dimensional wind field, considering the influence of all the unsteadiness that it brings. Both flight cases are done with an active control of proportional gains which were optimized for a certain wind scenario. Two methods for unsteady derivative estimation has been proposed and compared. Overall influence of unsteady effects on energy transfer between atmosphere and aircraft has been determined. Here we consider only longitudinal motion of the aircraft represented as a point-mass model. The derivation of the equations is followed by the previous works of Langelaan [START_REF] Langelaan | Biologically Inspired Flight Techniques for Small and Micro Unmanned Aerial Vehicles[END_REF] and Lawrence. [START_REF] Lawrance | Wind Energy Based Path Planning for a Small Gliding Unmanned Aerial Vehicle[END_REF] The influence of wind gradients on overall power of the aircraft in the function of climb angle have been done by Lawrence. [START_REF] Lawrance | Wind Energy Based Path Planning for a Small Gliding Unmanned Aerial Vehicle[END_REF] It is clear that aircraft have to perform certain manoeuvres according to knowledge of velocity field and control strategy in order to increase its power, absorbing the energy from the atmosphere. The kinematics of the aircraft is solved with Runge-Kutta ordinary differential equations solver 45 found in Matlab. The solutions of the system are velocity, climb and flight path angle, from which all other parameters of flight can be evaluated.

T -D -W sin γ = W g ( V + ẇx cos γ -ẇz sin γ) (2)
-L + W cos γ = W g (-V γ + ẇx sin γ + ẇz cos γ) (3)
M = θI yy (4)
The main objective of gust soaring is to increase the equation of power performing certain manoeuvres according to the current wind velocity perturbations. Energy gain can be achieved with optimal climb for negative component of vertical wind velocity which is known as up-draft and exploited by all glider pilots in thermals. Another way is through active control according to setting of wind derivatives.

Ėa

m = -gw z + qS m (C T -C D )V -V ( ẇx cos γ -ẇz sin γ) (5)
The control used in this work is based on set of proportional gains, which has been previously demonstrated by Langelaan. [START_REF] Langelaan | Biologically Inspired Flight Techniques for Small and Micro Unmanned Aerial Vehicles[END_REF] Simply, aircraft will react according to input information of the a priori known wind by flaps or horizontal tail. It should be pointed out that the gains have to be optimized with certain sensitivity to magnitude of the wind. Once obtained values for a certain wind profile will not be optimal for other with higher frequency and magnitude. This is due to the fact that every reaction has to outweigh drag increase in power equation.

δ = δ e δ f where δ e =      k 1 w z k 2 w x k 3 ẇx k 4 ẇz      and δ f =      k 11 w z k 12 w x k 13 ẇx k 14 ẇz      (6)
Wind velocity in two-dimensional model has been modelled as sinusoidal function. In this way we are able to recognize more obviously response of the aircraft and involve potentially beneficial unsteady effects depending on the frequency of posed wind profile.

w z = k sin(ωt) and ẇz = kω cos(ωt) (7)
The term endurance is related to a constant value very often calculated for cruise regime and represents ratio of available energy and power required. In the case of steady flight those parameters are quite obvious since velocity and altitude can be assumed. However in gusty flight the objective is to gain energy with time thus increasing specific power during the flight path. The benefit in endurance could be evaluated according to:

E nd E nd = t 0 1 - P P c dt (8)
The equation claims that by increasing the power with control P c (where P is power during the flight without energy harvesting control) we increase the magnitude of E nd , thus increasing the endurance of the aircraft.

Dynamic energy transfer between wind and aircraft

The goal is understanding of energy harvesting mechanisms through equations of motion and possible influence of some aerodynamic parameters. The flight of the airplane is considered to be symmetric and wind disturbances are coming from vertical plane. Once again the power equation can be transformed now revealing new aerodynamic parameter as shown in following equation:

Ėa m = -gw z -g 2n 3 ρ W S C D C 3/2 L -V ( ẇx cos γ -ẇz sin γ) (9)
The second term (representing power required) claims that the maximum energy transfer between wind and aircraft will be achieved at maximum value of the factor C 3/2 L /C D (corresponds to minimum sinking rate) for both stationary (n = 1) and non-stationary flight(n = 1). The term accounts as power required and should be minimized allowing higher levels of gained energy from wind power terms. It is also well known that in case of thermals (assuming stationary up-draft) glider pilot should estimate the velocity of the aircraft according to MacCready theory which refers to previous aerodynamic factor.

Further integration of wind power part reveals specific cases in which energy transfer can be maximized. It is obvious that presence of negative vertical wind component (up draft) increases the power of the aircraft. On the other hand considering wind derivatives in longitudinal plane we can differentiate some specific flight cases for maximization of energy transfer. Those specific flight cases for energy maximization have been obtained with detailed decomposition involving first and second order derivation of power equation parts related to wind.

Ėa,

ẇx m = -V ẇx cos γ = - ∂w x ∂x (V cos γ + w x )V cos γ - ∂w x ∂z (-V sin γ + w z )V cos γ (10) Ėa, ẇz m = V ẇz sin γ = ∂w z ∂x (V cos γ + w x )V sin γ + ∂w z ∂z (-V sin γ + w z )V sin γ (11)
Isolation of different wind parts and multiple derivation brings following extremes:

• Climbing/descending into positive/negative wind shear ∂wx ∂z would increase specific power of the aircraft where positive/negative horizontal wind shear would imply optimal angle γ = sin -1 wz±

√ 8V 2 +w 2 z 4V
• Optimal soaring strategy for power maximization in the case of positive/negative vertical wind shear ∂wz ∂x would be climb/descend for optimal climb/descend angle γ = cos -1 wx∓

√ 8V 2 +w 2 x 4V
• For a positive/negative linear vertical wind gradient optimal climb angle would be γ = sin -1 wz 2V /maximum climbing angle Note that the previous equations are developed with respect to the air-path system. On the other hand, the final form of power equation for inertial system is:

∂w x ∂z >0 <0 ∂w x ∂z >0 ∂w z ∂x >0 ∂w z ∂x
Ėi m = -q S m C D (V + w x cos γ -w z sin γ) -q S m C L (w x sin γ + w z cos γ) (12)
Wind gradients do not appear explicitly in the expression for total power expressed with respect to the inertial frame. However, time and spatially-varying winds will result in time varying forces (through changes in airspeed and angle of attack), so wind gradients will indirectly affect total power. In order to approach flight strategy in gusty environment with more realistic flight scenario, introducing three-dimensional effects, multi-point model has been developed. The model is consisted of three points representing two wind sensors on the wing and central mass-point. Since all points travel in a space grid of wind, trilinear interpolation 18 is applied taking into account all surrounding seed points of the mesh. This would provide alternative case of simulation to real flight test with equipped multihole probes capable of measuring three wind velocity components as has been demonstrated by Mohamed. [START_REF] Abdulrahim | Developement and Flight Testing of a Turbulence Mitigation System for Micro Air Vehicles[END_REF] Multi-point model allows estimation of unequal lift and drag distribution at each side of the wing, allowing instantaneous knowledge of necessary control action for coupling of roll and yaw moments. It also provides the information of gust length scale, frequency and magnitude on which can be based control strategy of the real flight, evading unnecessary actions in the case of small length scales. The equations of motion are now expanded taking into account lateral and roll motions. System of differential equations is once again integrated numerically using ODE45 solver with adapted time step to the wind velocity field.

T -D -W sin γ = W g ( V + ẇx cos ψ cos γ + ẇy sin ψ cos γ -ẇz sin γ) (13)
L -W cos γ = W g (-V γ + ẇx cos ψ sin γ + ẇy sin ψ sin γ + ẇz cos γ) (14)
L sin φ = W g (V ψ cos γ -ẇx sin ψ + ẇy cos ψ) (15)
The model assumes that the roll rate is directly controlled as an response from the unequal angle of attack along the wingspan. The two side points move according to the rotation of the aircraft angle ψ and emulate wing sensors as for example multi-hole probes. The information on wind velocity is then transformed into additional velocity and angle of attack seen by each side of the wing. Corresponding roll and yaw moment are then estimated with a lifting line theory demonstrated by Ringnes. 11 Ėa

m = -gw z + qS m (C T -C D )V -V ( ẇx cos ψ cos γ + ẇy sin ψ cos γ -ẇz sin γ) (16)
Control strategy has remained the same as for the longitudinal case. Once again the optimization of control gains has to be performed with the aim to increase power of the aircraft. Time histories of roll and yaw moment are showing the moment acting on wing that would be lost if no span-wise variation in turbulence were assumed. These moments are referred as added, and expressed as:

M roll,a = b/2 -b/2 l(y)ydy (17)
M yaw,a = b/2 -b/2 d i (y)ydy ( 18
)
where l(y) and d i (y) are sectional lift and induced drag as function of span-wise location y.

C. Unsteady aerodynamics

Generally, aircraft's lift and pitch moment coefficients can be represented by the MacLaurent series according to the following expressions:

C L = C L0 + C Lα α + C L α αl 2V + C L θ θl 2V + C L α αl 2 4V 2 + Control C Lδ δ + C L δ δl 2V + ... + ... (19)
C M = C M 0 + C M α α + C M α αl 2V + C M θ θl 2V + C M α αl 2 4V 2 + Control C M δ δ + C M δ δl 2V + ... + ... (20)
Aerodynamic derivatives can be divided in two groups. The first group can be be obtained by steady methods as for example modified vortex lattice method with included viscous effects while the other has to involve unsteady effects proposed in the following section. All of the previous authors omit the estimation of unsteady derivatives. Driven by the nature of turbulence, control activation is of high frequencies. Therefore, aircraft will perform rapid manoeuvres followed by the violent variation in angle of attack. Depending on the frequency and magnitude of pitching and heaving motion, trailing vortices are generated introducing unsteady behaviour of the aerodynamic forces. Theodorsen 21 modelling of aeroelastic airfoil has been published back in 1935, where aerodynamic forces are represented as a sum of non circulatory and circulatory effect for a thin section.

L = -ρb 2 (V π α + π ḧ -πbaα -V T 4 β -T 1 b β) -2πρV bC(k) V α + ḣ + b( 1 2 -a) α + 1 π T 10 V β + b 1 2π T 11 β (21)
Transfer function C(k) transforms a quasi-steady prediction of circulatory aerodynamic forces to unsteady values, introducing time lag effect. Since we are dealing with a time domain variations inverse Fourier transformation is applied for Wagner functions approximated by Jones. [START_REF] Jones | The Unsteady Lift of a Wing of Finite Aspect Ratio[END_REF] Those functions are also available for several aspect ratios, which enables comparison for various wing dimensions. The complete system of equations for non-linear response of a typical airfoil section exposed in the state-space formulation is done by Kholodar. [START_REF] Kholodar | The influence of a Nonzero Angle of Attack and Gust Loads on the Nonlinear Response of a Typical Airfoil Section with a Control Surface Freeplay[END_REF] Alternative way for estimation of unsteady derivatives is recently proposed by Gili. [START_REF] Gili | A New Approach for the Estimation of Longitudinal Damping Derivatives: CFD Validation on NACA 0012[END_REF] It involves unsteady Reynolds averaged Navier-Stokes equation solver coupled with rigid body motion available in Star-ccm+ software package. Attractive side of this method is that there is no need for re-meshing, deformable or interface mesh, which is highly cost effective in terms of time. The domain is the same as for steady computations with imposed motion law. Tested airfoild is SD2048 which belongs to the family of low-Re number foils. Structured hex mesh has been constructed with refined boundary layer modelling and number of elements have been chosen based on convergence criteria. The time step chosen for unsteady simulations is of the order 0.5E-4. Since Theodorsen model is developed for flat plate, certain modifications had to be done for adequate comparison with airfoil. In order to fit curves of resulting lift coefficient, translation for C Lo and rotation for C Lα obtained from steady simulation had to be implemented for theoretical model of Theodorsen. After implementation of corrections certain gap in lift slope was found due to the viscous effects modelled by CFD. Nevertheless, acceptable comparison has been found, both for harmonic pitch and coupled motion of pitch and heave. Once the unsteady coefficients are determined we proceed to lift representation in energy harvesting simulations according to the frequency of the gust. The previous effort was realized in order to model realistic case of lift evolution during energy transfer between atmosphere and aircraft followed by aggressive manoeuvres. For rapid manoeuvres, it is important to include added-mass terms to account for the reaction force due to the mass of the fluid that is accelerated by the airfoil. Additionally, one must include induced circulation around the airfoil due to wake vorticity. Moreover, importance of taking into account unsteady effects has been justified with improved harvesting cycles. Implemented unsteady derivative coefficients featuring time-lag result in increased lift coefficient (as shown in Fig. 6) during harvesting cycles. Following large accelerations lift coefficient is strongly influenced by added-mass forces according to Brunton. [START_REF] Brunton | Empirical state-space representation for Theodorsen's lift model[END_REF] Those increments can bring additional positive work on aircraft, in comparison to steady model, increasing harvested energy. Also it should be pointed out that frequency response is qualitatively in relation to the position of rotating point. If the location of rotation point goes afterwards of mid-cord the effect of added-mass with higher frequencies becomes negative. With a positive steps in angle of attack lift would move in negative direction due to negative added-mass value outperforming circulatory forces.

V. Structural dynamics

T his chapter describes the physics behind passive method of energy extraction by means of energy transfer from gust without autopilot interference. The presented study aims at discovering the aero-elastic effect of wing on overall energy gain in energy extraction flight regime. Related work has been performed by Mai 28 on energy reteieval for sailplanes penetrating a gust and Ironside 14 on wing drag reduction due to structural dynamics. Strategy also brings an insight for structural tailoring during the process of UAV design for passive control. The goal is to determine the impact of wing motions on flight dynamics and potentially quantify the benefits on energy gain. Two coupled motions are considered in function of time provoked by variable loading. Firstly, flapping due to variable lift force and secondly, torsion due to the distance between aerodynamic and shear center.

A. Elastic wing

The wing was divided into arbitrary number of spanwise elements. Each section is represented by its own cord, center of shear, center of gravity and aerodynamic center. Each wing element is assumed to have a given mass acting on the local center of gravity and surface. Aerodynamic and mass forces are assumed to act on aerodynamic center and center of gravity, respectively. Twisting was assumed to occur about the elastic axis. In structural dynamics, wings are usually represented as thin beams [START_REF] Bishop | On coupled bending and torsional vibration of uniform beams[END_REF] and therefore are represented in following state-space form. Obviously the size of the following system of equations of motion depends on the number of elements chosen for wing. Composition of the matrices in function of chosen number of nodes was described by Panzer. 26

F 0 0 M Ẏ Ÿ (t) = 0 F -K -D Y Ẏ (t) + 0 L (t) (22)
Where F , M , K, D and L are identity, mass, stiffness, damping and load matrices, respectively. The local angle of attack induced by a wing motion acting on a single section is computed with following expression:

α i,str = θ xi - żi V + δ ci θ xi V (23)
The previous value of induced angle of attack is used as an input of aerodynamic model which calculates the lift distribution due to the geometry of the wing itself. Wing model was taken from UAV DT-18 shown in Appendix. The equations of motion for wing and aircraft flight mechanics were solved using Runge-Kutta scheme and input of adaptable mass and stiffness matrix. Sectional lift and moment equations due to structural dynamics were computed with following formula:

M istr = q o S i aδ ei θ xi - żi V + δ ci θ xi V and L istr = q o S i a θ xi - żi V + δ ci θ xi V ( 24
)
where δ ei is the distance between the aerodynamic center and center of shear and δ ci is the distance between center of gravity and downwash control point. It can be concluded that main parameters that contribute to wing deformations are chosen isotropic material represented by Young Elastic modulus E, length of the wing, spar moment of inertia and sectional parameters as aerodynamic, shear and gravity center. Despite the fact n i (t+dt)

n i (t) γ i+1 SC i AC i θ i C i C i SC i γ i β i+1 β i θ i θ i+1 w i w i+1 v i v i+1 u i u i u i+1 n i i i+1 CG i Figure 7. Deformation principle.
that wing deformation are of small scale and that their contribution is of low order on energy extraction, passive approach still represents very attractive strategy, since it does not require any knowledge of the wind field.

VI. Results

T wo small UAVs have been tested with aerodynamic properties exposed in Appendix. The goal was to quantify the level of energy harvesting and recognize aerodynamic parameters which could enhance flight strategy.

A. Longitudinal simulation results

Simulations in vertical plane have been performed with imposed sinusoidal wind field characterized with frequency and amplitude shown in Fig. 11. Two flight cases have been examined from which, one with active control of proportional gains and the other without control. Unsteady derivatives have been chosen according to the frequency of imposed wind field. The values of proportional gains have been optimized for one frequency and magnitude of presented flight case. The implementation of correct unsteady derivatives showed to be with a slight positive impact on energy harvesting. Also, it should be pointed out that, the gains optimized ones are not any more optimal for different frequency or magnitude, which implies different control strategy in case of real flight test.

Two different aircraft served as aerodynamic model for simulations. On one hand XC glider of twice longer wing-span, higher lift-to-drag ratio and higher mass, while on the other side DT-18 of higher control efficiency (see Appendix for aircraft properties). The results show evident energy gain (see Fig. 11) for both vehicles. It is interesting to notice that XC glider being more aerodynamically efficient, since of higher lift-to-drag ratio, achieves higher level of energy after 100s of flight. However, in case of DT-18 flight, energy harvesting is on higher magnitude for about 15%. One of the missions of this paper was to detect which parameters would contribute to such a divergence in energy gain. One of the parameters recognized for strong influence on energy harvesting was certainly lift slope coefficient (C Lα ). It was found out that 15% higher value of lift slope brings 12% more energy after 100s and opposite for lower values. Another parameter affecting is control effectiveness represented by the coefficient (C Lδe ). More efficient elevator of 25% in comparison to original plane would bring 8% more energy after 100s of flight. Moreover, flying at higher zero angle of attack lift coefficient (C Lo ) brings significant benefit in energy gain. For 10% higher value of C Lo the aircraft would achieve 20% more energy with active control.

Optimal definition of previously mentioned parameters should be objective in the process of UAV design for energy extraction and long endurance flight. They indicate how important is effective gain in lift during manoeuvrable flight within the energy harvesting cycles. As usual drag reduction is highly beneficial as it directly affects power equation. Gaining energy with gust soaring is a constant battle between beneficial part of wind components and drag increment due to reaction as described by power equation. Topology of final energy achieved in function of control gains has been built in respect to constraints of the motion during energy extraction. The optimization problem has been set with objective function of maximizing final energy achieved using Non sorting genetic algorithm II (NSGA2) and Sequential least squares algorithm (SLSQP) based on quasi-Newton method. Control strategy in presence of only vertical wind was based on two control gains K 1 and K 2 multiplying a priori knowledge of w z and ẇz respectively. The case with horizontal profile would involve 2 more coefficients, respectively. Despite the fact that genetic algorithm brings highest energy harvesting result, computation time is relatively high. This was important due to the fact that proportional control has to be optimized for each frequency and magnitude of wind profile. On the other hand, attractive algorithm based on gradient method achieves faster convergence, although strongly dependent on initial values given to optimizer. With this purpose rough topology has been constructed shown in Fig. 8 in order to have more precise idea about initial values. The difference between two optimization methods is negligible and goes into favour of SLSQP method with proper initiation and Jacobian matrix definition, thus reducing significantly computational time. Both optimization algorithms included are used in openMdao package for Python. With the aim of more realistic flight case, three dimensional flight dynamic simulations have been performed. Flight environment was three-dimensional wind field of defined spatial and temporal resolution. The vehicle was modelled as multi-point system, now collecting additional information about wind velocity from side sensors located on the wing. Based on assumption that a single point is adequate for angle of attack estimation of whole wing side, roll and yaw moments are anticipated. Energy harvesting mechanism once again proved to bring significant benefits, even in such a complex wind field, bringing around 20% more energy after 60s of flight with active control. Trajectory and energy results shown in Fig. 12 confirmed behaviour of longitudinal flight. Once again XC plane proved to be more efficient in terms of aerodynamic performance, while on the other hand DT-18 more efficient in energy harvesting. Side points which represent wing sensors (ex. of multi-hole probes or pressure based holes on the wing used for the real flight test) are moving according to wing kinematics. The information provided by those points (see Fig. 13) are then transformed into roll and yaw moment provoked by the unequal wind velocity seen by each side of the wing. In the real flight test campaign those information would be used for the activation of ailerons, coping unwanted actions of the plane. They would also serve as a decisive mechanism triggering reaction for energy harvesting due to possible unequal length scales of gust.

As artificial flying environment depends on his creator, thus coherence level is subject of definition in tree-dimensional wind field used in this work. It is defined with a specific function that can be found in wind section of IEC-61400 standard.

C oh (r, f ) = e -12((f r/V hub ) 2 +(0.12r/L C ) 2 ) 0.5 (25)
The chosen parameters result in higher coherence of horizontal velocity (blue line Fig. 9 (b)) in comparison to vertical one. This is the reason why we are coping more roll than yaw moment in our simulations (shown in Fig. 13) as the lift increment is higher due to the deviations in angle of attack than in horizontal velocity.

C. Elastic wing in flight

The program solves the equations of motion, both for the wing and aircraft this way taking into account structural dynamics. The simulations have been done for the sinusoidal vertical wind field. Sectional force and moment are computed according to the local angle of attack for each wing section, provoked by the aircraft manoeuvring, wind itself and structural dynamics. As explained in Chapter V the angle of attack induced by wing motion is due to velocity of vertical deformations, local torsion and its velocity, coupled all together. The following figure shows exaggerated both torsion and coupled deformations due to the inequality in axis scales. Although, the model is adapted for arbitrary choice of wing sections, the recommended number would be 10 for the existing wing model, concerning the computational time.

The results show beneficial effect of structural dynamics on overall energy retrieval as shown in Fig. 10 (c). It is interesting to point out that structural dynamics of bending contributes more than torsional on final energy gain. Due to the previous statement, aero-elastic tailoring should be subjected in the process of wing sizing as it contributes to energy gain. Moreover, the parameter δ ei representing the distance between aerodynamic center and center of shear contributes to energy harvesting with its growth. As center of shear depends on sectional geometry which is defined in aerodynamic design, possible manipulation could be achieved with position of different elements inside section as for example spar. Wing elasticity in service of energy harvesting is particularly interesting technique as it does not require any knowledge of the wind field. However, its contribution is much lower in comparison to active control.

VII. Conclusion

T his paper demonstrates the feasibility of gust soaring for a small UAV, theoretically quantifying the gain in energy through numerical simulations. Beside the unrealistic sinusoidal wind profile, although easy to manipulate, complex and stochastic three-dimensional wind field has been used as a flight environment where realistic wind profiles were used as a reference for artificial three-dimensional wind field generation. Model of unsteady aerodynamics has been implemented in order to simulate effects provoked by aggressive and rapid manoeuvres compliant to stochastic nature of wind. Three-dimensional effects on the wing due to unequal wind experienced along wing-span has been transformed into roll and yaw moment for the multipoint model. Finally, passive approach was examined, particularly interesting as it does not require any knowledge of the wind, showing benefits and revealing the significance of careful aero-elastic tailoring with emphasis on energy extraction.

Inevitable future contribution on the subject of gust energy harvesting would be definition of the design process according to relation between the aircraft geometry and energy gain. Particular interest will be to define an optimal design strategy for a UAV which would more effectively extract energy from the gust.

Gust soaring theoretically proved to be very promising flight technique for endurance enhancement. However, in order to reach full understanding about refined kinematics and energy exchanges, it seems necessary to perform the flight test campaign including wind measurement system equipped on the aircraft. It would bring requisite knowledge of the real wind length scales and challenging control strategy.

Figure 1 .

 1 Figure 1. Three-dimensional wind field.

Figure 2 .

 2 Figure 2. 2D flight physics off mass point model.

Figure 3 .

 3 Figure 3. Flight cases for power maximization in wind shear (Note that for most common ratio between V, wx and wz angle γ is close to the value of ±45 • )

Figure 4 .

 4 Figure 4. 3D Flight physics of multi-point mass model.

Figure 5 .

 5 Figure 5. Unsteady derivative estimation comparison between Theodorsen and URANS simulations.

Figure 6 .

 6 Figure 6. Unsteady vs steady cycle during sinusoidal vertical gust

Figure 8 .

 8 Figure 8. Optimization topology for sinusoidal gust

Figure 9 .

 9 Figure 9. Tree-dimensional flight

Figure 10 .

 10 Figure 10. Wing deformation within vertical sinusoidal wind gust (wz = 2 sin(2t)).

Figure 13 .

 13 Figure 13. Roll and Yaw moment estimation for XC glider

Figure 15 .

 15 Figure 15. DT-18

Table 1 .

 1 Optimization constraints

	Parameters of control gains optimization
	Optimization objective: Final energy achieved E a,f inal	Maximize
	Velocity V range	[10 : 30 m/s]
	Initial V	18 m/s
	AoA α no stall range	[-6 : 12 • ]
	Flight path angle θ range	[-50 : 50 • ]
	Optimization package OpenMdao with 100 generations for NSGA2

Table 3 .

 3 DT-18 glider by DelairTechC 1 α 5 + C 2 α 4 + C 3 α 3 + C 4 α 2 + C 5 α + C 6

		Aerodynamic properties of DT-18
	Parameter	Value	Description
	M	1.9 kg	Mass of the aircraft
	b	1.8 m	Wingspan
	l	0.18 m	Mean aerodynamic cord
	S	0.248 m 2	Surface of wings
	I xx	0.184 kgm 2	Roll moment of inertia
	I yy	0.12 kgm 2	Pitch moment of inertia
	I zz	0.66 kgm 2	Yaw moment of inertia
	λ	13	Wing aspect ratio
	(C L /C D ) max 15	Lift to drag ratio
	C Lo	0.4	Lift coefficient for α = 0
	C Lα	6.37 /rad	Lift slope
	C Lq	13.9254 s/rad
	C Lδe	0.6188 /rad
	C Lδ f	2.158 /rad
	C D		Drag coefficent
	C Dδe	0.02292 /rad
	C Dδ f	0.0286 /rad
	C M q	-24.68 s/rad
	C M α	-3.22 s/rad
	C M δe	-2.4977 /rad
	C M δ f	-0.3953 /rad

Downloaded by UNIVERSITY OF COLORADO on January 14, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2017-1630

Appendix

A. Aircraft properties XC glider

The aircraft model used in simulations is XC glider designed by Mark Drela. There are numerous ways for estimation of steady aerodynamic coefficients, both analytical or numerical. The coefficients presented in Table 2 are obtained with several software packages based on vortex lattice method. A new software which enables implementation of viscous corrected polars for section airfoils generated from XFOIL is called openVogel. For the sake of comparison and curiosity another software is used called Tornado.