Dynamic nonlinear aeroelastic model of a kite for power ... - AWESCO [PDF]

Dynamic nonlinear aeroelastic model of a kite for power generation Allert Bosch,∗ Roland Schmehl,† Paolo Tiso‡ and Daniel Rixen§ Delft University of Technology, Delft, The Netherlands.

This paper presents a numerical modeling approach for the flight dynamics and deformation of tethered inflatable wings which are central components of many contemporary airborne wind energy systems. A geometrically nonlinear finite element framework is used to describe the large quasi-static deformations caused by changing bridle line geometry and by varying aerodynamic loading. The effect of the external flow is described in terms of discrete pressure distributions for the different wing sections. The empirical model takes into account the shape parameters chord length, camber and thickness per section and is derived by fitting precomputed data from computational fluid dynamic analysis. To reduce computation times, local dynamic deformation phenomena are neglected. Each integration time step, the steady aerodynamic loading is determined first and then used to update the static equilibrium shape of the wing. This static aeroelastic model is embedded in a dynamic system model which includes the tether, bridle lines and kite control unit. The iterative approach can accurately describe bending and torsion of the wing, which contribute to the aerodynamic steering moments. The presented approach is complemented with a flight controller and used to simulate figure-eight flight maneuvers of a leading edge inflatable tube kite used for traction power generation.

I.

Introduction

Airborne wind energy (AWE) systems use tethered flying devices that can operate at higher altitudes than conventional wind turbines. A main advantage is the expected significant increase in capacity factor, the ratio of actually generated to potentially feasible power output, which is a result of the more constant and stronger wind [1–3]. Next to this, replacing the rigid tower of a wind turbine by a lightweight tensile structure translates directly into lower investment costs and a lower environmental footprint. The reduced visual and acoustic impact is an advantage for installations in ecologically sensitive areas, while the low weight and compact dimensions are particularly suitable for mobile deployment. Among the various concepts that have been devised and tested over the past decade, a particular focus has been the conversion of the traction power of a tethered inflatable wing into electricity [4–7]. As depicted in Fig. 1, such a kite power system is operated in pumping cycles consisting of two alternating phases. During the reel-out phase the kite is flying

Wind

Reel-out phase: energy generation

Reel-in phase: energy consumption

Figure 1. The two phases of a pumping cycle [7]

fast crosswind maneuvers, e.g. figure-eights or circles, to create a high ∗ Researcher, ASSET Institute, Kluyverweg 1, 2629HS, Delft, The Netherlands, [email protected]. † Associate Professor, ASSET Institute, Kluyverweg 1, 2629HS, Delft, The Netherlands, [email protected]. ‡ Assistant Professor, Dept. Precision and Microsystems Engineering, Mekelweg 2, 2628CD, Delft, The Netherlands, [email protected]. § Professor, Dept. Precision and Microsystems Engineering, Mekelweg 2, 2628CD, Delft, The Netherlands, [email protected].

1

traction force. This traction force is driving an electrical generator on the ground. During the reel-in phase the angle of attack of the kite is decreased in order to reduce the traction force. As consequence, the kite can be retracted to a lower altitude for only a fraction of the energy that was generated in the reel-out phase. The kite is steered and de-powered by means of a kite control unit (KCU) which is mechanically connected to the wing by a bridle line system and to the ground station by a tether. A wireless link is used for remote-control and streaming of airborne sensor data. In contrast to the rotor blades of a conventional wind turbine, an inflatable membrane wing is lightweight and highly flexible. Its shape is determined by the interaction of aerodynamic wing loading and pressure distribution in the tubes on the one hand and internal forces in the membrane fabric and bridle lines on the other hand. Due to its low inertia, the wing responds generally very sensitively to control commands, fluctuations of the apparent wind velocity and other changes of geometry or boundary conditions. This is one of the reasons why robust automated flight control is a challenging task [8]. In essence, the strongly coupled aerodynamics and structural dynamics of the wing constitute a challenging fluid-structure interaction (FSI) problem which interacts with the flight dynamics. On the other hand, the boundary conditions of the FSI problem continuously change as the relative flow conditions vary along the flight trajectory. Computational analysis has confirmed that the spanwise torsion of the wing plays an important role in the steering behavior of the wing [9]. Other aeroelastic phenomena such as spanwise bending oscillations or wing tip flapping can be induced by the interaction of aerodynamic and structural forces. Many studies with focus on system performance or flight control consider only point mass or rigid body flight dynamics and do not account for effects of wing flexibility [10–13]. This significantly limits the accuracy and prediction capability of the simulation model because the physical mechanism of steering by means of line displacement can not be represented adequately. Standard approaches for aircraft flight control, such as system identification, are difficult to use because the deformation of the wing strongly depends on the time history of the flight maneuvers [14]. One possible improvement is to include an empirical correlation between steering input and induced flight dynamic response of the wing, such as the correlation between steering line displacement and yaw rate proposed in [15–18]. However, despite its mathematical simplicity, such a black box model requires a statistical evaluation of a sufficient volume of flight data. Design changes generally require a renewed statistical evaluation and, accordingly, such a model has limited practical value in the early stages of development when prototype flight data is often not available. An attempt to include mechanistic model components is proposed by [19], assuming a correlation between steering line displacement and induced spanwise wing torsion. This complex three-dimensional deformation mode is idealized as a shear distortion of the projected wing

2

geometry, described by a single deformation parameter, the wing tip displacement, which depends linearly on steering line displacement. The focus of this study is on the coupling with a panel method to determine the aerodynamic surface pressure distribution on the deformed wing. The central assumption is thus potential flow, however, most relevant flight situations are at a high angle of attack, which implies separated flow, which is an inherently viscous phenomena. The practical use of this method for quantitative analysis is limited. A similar approach is presented in [20]. To systematically improve the wing design and to optimize its operation under a wide range of wind conditions a more detailed understanding of the flight dynamic behavior is indispensable. This is a major motivation for the development of mechanistic models that incorporate deformation on the basis of structural degrees of freedom (DOF). In order to accelerate simulation times, small-scale FSI phenomena, such as canopy flutter which is caused by local pressure fluctuations and the low inertia and stiffness of the wing, are typically neglected as the impact on the overall flight dynamics is minor. Macro-scale FSI phenomena, such as wing torsion and bending, have a major effect on the flight dynamics and are taken into account. Some studies model the aeroelasticity of the wing on the basis of a particle system representation [21], a multi-plate representation [22] or a multibody representation [9, 23]. However, these type of mechanistic models generally imply high degrees of idealization, especially in defining the mechanical interaction between the discrete structural elements. Accordingly, the required empirical parameters vary with the specific wing design and consequently need to be determined experimentally on a case-by-case basis. Closer to the physical reality is a consistent finite element (FE) representation of the membrane wing. In a recent study, several 10000 triangular membrane elements are used to discretize the fabric material of a leading edge inflatable (LEI) tube kite together with the line elements for the bridle system [24]. The internal aerostatic loading of the pressurized tube chambers is described by a uniform pressure model. An exterior aerodynamic loading of the wing is not considered. Such level of physical detail allows for predictive simulation of the wing deformation requiring only the material characteristics of the membrane fabric and line materials. However, the high-resolution discretization inevitably resolves dynamic phenomena with short time-scales, which reduces the numerical integration time step and thus increases the computational time. The static deformation of a ram air parachute in flight is analyzed in [25]. The method iteratively couples a FE representation of the elastic membrane wing with a vortex lattice method (VLM) to determine the surface pressure distribution generated by the exterior flow. An improved approach is presented in [26], coupling a FE representation of the ram air wing with a panel method to simulate the exterior aerodynamics. The inviscid flow solver includes empirical corrections for several wing components and a simple boundary layer model to account for flow separation. From this review of existing research it is apparent that an accurate and efficient mechanistic model of a flexible tethered wing would be a valuable contribution to the present AWE development activities. The objective of this paper is to derive such a model, covering the flight dynamics of the wing as well as its quasi-static deformation under aerodynamic load. The mechanistic approach provides detailed insight into the behavior of the wing during specific flight maneuvers and establishes a generic link between steering line displacements and large-scale dynamic system response. Small-scale phenomena such as local fluttering of the canopy are not considered. This approximation drastically reduces the overall computational effort. A basic flight controller is included in order to simulate figure-eight flight maneuvers that are used for many AWE Systems. The presented approach is discussed further in [27, 28]. Following this introduction, Sect. II describes the computational approach and the individual model components, followed by a description of the flight controller in Sect. III. The simulation results for figure-eight flight maneuvers are discussed in Sect. IV and the conclusions are drawn in Sect. V.

II.

Computational approach

A common approach for the numerical simulation of complex FSI problems is the alternating solution of fluid and structural dynamic problem parts within an iterative coupling procedure. Since larger deformations can significantly affect the fluid flow and the resulting surface pressure distribution two-way coupling is of crucial importance. This can be achieved by exchanging updated structural displacements and fluid pressure at the interface during the coupling iterations. To simulate the deformation of the tethered wing along a crosswind flight trajectory with continuously changing relative flow conditions and bridle line geometry the FSI problem is embedded in a dynamic system model. Since the objective of the present study is the description of the quasistatic large-scale deformation response, a static FSI problem is solved by coupling a static structural model with a steady aerodynamic load model. Small-scale deformation phenomena with short time scales (high frequencies) are effectively filtered out by this approach. As mentioned in the introduction, these phenomena generally have only a marginal effect on the dynamic behavior of the wing under nominal flight operation. Figure 2 illustrates the iterative procedure comprising the three distinct component models. The static structural model is based on a Dynamic system model

vw

zw yw xw ∆t

Displacements

Forces

Static structural model zk yk

xk Displacements

Surface load

Steady aerodynamic load model Static FSI problem

Figure 2. Iterative coupling procedure.

finite element discretization of the wing. The canopy is represented by triangular shell elements, whereas the pressurized tubular frame is represented by conventional beam elements. Compared to other studies [24], the discretization is very coarse such that small-scale phenomena are not resolved. Wrinkling, for example, has an important influence when a beam element is bent strongly and eventually buckles. However, such severe off-design conditions are excluded from the scope of the present analysis. The steady aerodynamic load model provides a discrete surface pressure distribution for each wing section as a function of its local relative flow conditions and shape deformation. The approach is adapted from [9, 29] and matches the level of detail and computational effort of the structural model. The dynamic behavior of the complete system consisting of wing,

3

bridle lines, KCU and tether is described by a particle system model. Bridle lines are represented as linear spring-damper elements, whereas the tether is described as a distance constraint. The inertia of the system components is represented by point masses at the KCU and at the bridle line attachments on the wing. For each dynamic integration step, the structural model provides the resultant wing forces at the bridle line attachments. These forces, which constitute aerodynamic and internal structural components, are applied as external loads to the particle system model. Following the integration step, the updated positions of the point masses are returned to the structural model and treated as displacement boundary conditions. Each of the component models depicted in Fig. 2 uses an individual reference frame. The dynamic equations of the particle system are formulated in the wind reference frame (xw , yw , zw ), an inertial frame which has the xw -axis aligned with the wind velocity vector vw and the zw -axis perpendicular to the ground surface. The structural deformation of the wing is analyzed in a reference frame (xk , yk , zk ) moving with the kite. The flight velocity vector vk and the angular velocity vector ωk describe the translation and rotation of this kite-fixed reference frame with respect to the wind reference frame. Together with the wind velocity vector vw the vectors vk and ωk constitute the state vector X of the wing. The aerodynamic loading is determined in reference frames (xki , yik , zik ) that are moving with the respective wing sections. A.

Static structural wing model

The deformation of the wing is described by the displacements of the   structural nodes. The nodal displacement vector q = qf , qb comprises two distinct partitions: qf , containing the degrees of freedom and qb , containing the prescribed displacements of the bridle attachment points. For a specific flight condition and bridle line geometry the static deformed shape of the wing at an arbitrary time is determined by the equilibrium of internal elastic forces Fi and external aerodynamic forces Fa , which is described by a system of nonlinear equations,  f  F (q) = Ffa (q, X),     ib (1) Fi (q) = Fba (q, X) + Fb ,     q0 , qb , X prescribed, where q0 are the known displacements at the previous integration step and Fb are the unknown reaction forces at the nodes with prescribed displacements. The superscripts f and b denote the components of the nodal force vector at the free and prescribed displacement nodes. The kite used in the present study is shown in Fig. 3. The displace3.20 m

1.93 m

Leading edge Strut

2.94 m

Canopy Right tip A Tip beam

zk yk

B xk

rk

5.51 m 5.80 m Steering line Power line

C D Left tip 0.634 m

Figure 3. FE mesh, components and nominal dimensions of the North Rhino 16 m2 kite.

ments are prescribed at the four bridle line attachment points, which are

the corners A, B, C and D of the wing. These four points are also used to define the reference frame (xk , yk , zk ) for the structural analysis. The origin rk of this frame is positioned on the geometric midpoint between the two leading edge attachment points. The xk -axis is pointing to the midpoint between the two trailing edge attachment points. The yk -axis is pointing towards the leading edge attachment point on the right wing tip. The system of equations (1) is nonlinear in the displacements q, which has three different reasons. Firstly, large displacements cause a coupling between out-of-plane bending and in-plane stretching. This is accounted for in the adopted von K´arm´an kinematic description, which establishes a quadratic relation between the strains and the displacements, significantly changing the stiffness of the system under operational loads. The von K´arm´an kinematic model is adequate for the expected degree of deformation in typical flight situations. However, the described procedure is general and can include more complex kinematic models. The stiffening effect is essentially due to the reorientation of the canopy membrane when subjected to transverse pressure. The reorientation causes tensile stresses that are resisted by the much stiffer in-plane rigidity of the canopy Secondly, the bending of inflated tubes affects the tube cross section which in turn changes the bending stiffness. At larger deflections, this nonlinear geometrical effect can cause buckling of the tube, which essentially is a sudden decrease of bending stiffness. In order to model this phenomenon accurately, a very fine structural mesh would be required to capture the associated local wrinkling of the membrane. Since the objective of the present study is computational efficiency, the inflated tubular components are discretized by conventional two-node beam elements, i.e. omitting a discretization of the tube cross section. The nonlinearity is incorporated by deflection-dependent cross-sectional properties that are calibrated with experimental data [9]. Thirdly, the aerodynamic load distribution Fa (q, X) depends nonlinearly on the shape deformation of the wing, as explained in the following section. Severe nonlinearities are introduced by flow phenomena such as flow separation which, in the extreme, can lead to stalling of the wing. Due to the inherent nonlinearities, the system of equations (1) has to be solved iteratively. A Newton-Raphson scheme is adopted to solve the first block of equations in system (1) for the unknown nodal displacements qf while keeping the aerodynamic loads Fa constant. After convergence of the structural loop the aerodynamic loads are updated for the new shape of the wing and the displacements are recalculated. This procedure is repeated until the prescribed convergence criteria on displacements and force residuals are met. The reaction forces Fb occurring in the bridle line attachment points are then calculated by means of the second block of equations in system (1). The specific wing analyzed in the present study is a North Rhino kite with a total surface area of 16 m2 . The FE discretization of this kite is depicted in Fig. 3, including a description of essential components and global dimensions. This particular kite was chosen since it represents a typical kite used for AWE systems. The C-shaped canopy is supported by a pressurized tubular frame designed as a leading edge tube with five connected strut tubes. Each of the wing tips incorporates a thin stiff beam with a front and rear attachment of the bridle lines. The two front bridle lines, attached at points A and C, are denoted as power lines because they transfer most of the traction force generated by the kite to the tether. The two rear bridle lines, attached at points B and D, are denoted as steering lines because they are actively controlled to steer, power and de-power the kite. The trailing edge incorporates a thin wire to prevent over-stretching and flapping of the canopy. The canopy is made from ripstop fabric material and attached to the top of the tubular frame. The dimensions specified in Fig. 3 refer to the unloaded and undeformed state of the wing, with nominal angles of attack of αl = αr = 6.5°and αm = 12.7°at the tips and the midwing sections, respectively. The canopy is discretized by nonlinear three-node triangular flat shell elements combining a membrane element for the in-plane forces as proposed by [30, 31], and a bending element for the out-of-plane forces as proposed by [32], amounting in a total of 18 DOF. The ripstop fabric material is modeled as isotropic linear material with a Young’s modulus E = 1250 MPa, a Poisson ratio ν = 0.3 and a thickness h = 0.08 × 10−3 m. The positioning of the finite element nodes on the wing takes into account the method of determining the discrete

4

aerodynamic loading discussed in the following section. These forces are based on a discretization of the lifting surface in spanwise and chordwise directions which accordingly leads to the triangle strip mesh structure depicted in Fig. 3. This mesh is sufficiently fine to resolve the large-scale deformations of the wing while keeping the number of DOF at a reasonable level. The shorter chord in the tip areas leads to a finer mesh which is favorable because these areas experience the largest deformations and have the largest impact on the steering behavior. The pressurized leading edge and strut tubes, the trailing edge wire and the tip beam are discretized by nonlinear beam elements. This type of element is based on classic Bernoulli beam theory without shear deformation and has total of 12 DOF [33, 34]. The geometrical nonlinearity related to the varying cross section of inflated tubes under bending load is taken into account by means of a bending stiffness EI which, for a certain choice of tube material, depends on the inflation pressure p, the diameter d and the tip deflection q of the tube. The correlations EI(p, d, q) have been derived in [9] by fitting to experimental data. In total, the FE model comprises 107 beam elements, 360 shell elements, 222 nodes and 1332 DOF. B.

Aerodynamic wing loading model

To determine the aerodynamic surface loading the wing is partitioned into spanwise sections i = 1, . . . , n which are represented by twodimensional membrane airfoils with load points j = 1, . . . , m. As illustrated in Fig. 4, the nodes of this aerodynamic load mesh coincide with the nodes of the structural surface mesh. The approach adopted Fi2 L

Fi1 L

Fi3 L

Fi4 L

via,xz zik

which indicates the contributions due to translation and rotation of the kite-fixed reference frame (xk , yk , zk ) and due to deformation of the wing structure. The deformation velocity ui describes the velocity of the frame (xki , yik , zik ) relative to the frame (xk , yk , zk ) which can be decisive for providing aerodynamic damping. Since the quasi-static deformation model does not provide nodal displacement velocities, ui is approximated as a finite difference using the displacements of the leading and rear edge nodes of each wing section " # qim − qim 1 qi1 − qi1 0 0 ui = + . (3) 2 ∆t ∆t This average deformation velocity is well suited to describe the contribution of large-scale deformations which cause a translation or rotation of the wing section relative to the kite-fixed reference frame and which thus affect the local relative flow conditions. The spanwise component of via has a negligible effect on the aerodynamic loading. Consequently, only the projection into the xz-plane of the wing section via,xz = via − (via · eiy )eiy

c= b

d

i1 + r − rim

2

(

c

κ = max

j=1,...,m

Wing section

1 ρ CL c w (via,xz )2 , 2 1 FDi = ρ CD c w (via,xz )2 , 2 1 i M = ρ CM c w (via,xz )2 . 2 FLi =

w

Figure 4. Wing section i in side view (top) and perspective view (bottom).

from [9, 29] has the advantage that the fluid-structure coupling does not require computationally expensive interpolation of displacements and forces at the interface. The aerodynamic surface loading is determined for each wing section individually, by means of computational fluid dynamics (CFD), taking into account the relative flow conditions and the instantaneous shape of the specific section. Since the flow time scales are much shorter than the large-scale deformations of the wing, only the steady aerodynamic loading is of interest. Since the flow around a membrane airfoil with leading edge tube is characterized by significant flow separation, even for low angle of attack, the two-dimensional flow simulations are based on Reynolds-Averaged Navier-Stokes (RANS) equations using a k − ω turbulence model. Because such analysis is computationally extremely demanding it is not performed as part of the numerical integration procedure, but used to precompute a correlation framework. This approach allows an efficient evaluation of the surface loading for given relative flow conditions and shape of the airfoil. The local relative flow conditions are described in reference frames that are moving with the respective wing sections. As indicated in Fig. 4, the origin of wing section reference frame (xki , yik , zik ) is positioned at the leading edge node, with the xki -axis pointing towards the rear edge node and the yik -axis defined perpendicular to the airfoil. Defining the node

) eiz · (ri j − ri1 ) . c

(6)

The correlation framework is based on the following representation of the integral aerodynamic load on wing section i

xki

Finite element

(5)

and the camber as

α

yik

(4)

is used to characterize the local relative flow conditions. The angle of attack α is measured between vector via,xz and the the airfoil chord line. The airfoil geometry is described by the chord length c, the camber κ = b/c and the thickness t = d/c. In the discrete model, the chord length is evaluated as

Fi5 L Fi6 L

d

position vectors ri j as relative to rk , the local apparent wind velocity via is evaluated as via = vk + ωk × ri1 + ui − vw , (2)

(7) (8) (9)

The aerodynamic coefficients CL , CD and CM are regarded as correlations of the non-dimensional parameters κ, t and α. For all discrete parameter combinations within the range of possible operating conditions CFD simulations of the relative flow problems are performed to calculate the resulting aerodynamic coefficients. The fitted continuous functions are valid for −20° < α < 20°. Outside of this range, CL and CD are described by the empirical functions proposed in [35]. During the numerical integration procedure, the integral load characteristic, represented by a combination of aerodynamic coefficients CL , CD and CM , has to be expanded back to a consistent set of aerodynamic nodal forces Fiaj = FiLj + FiDj . In order to do so, the force components FiLj and FiDj are distributed over the airfoil nodes j = 1, . . . , m in such a way that the prescribed airfoil moment M i is matched. The employed varying weighted approach is explained in detail in [9, 23]. C.

Dynamic system model

In order to simulate complete flight maneuvers the static aeroelastic model of the wing is integrated with a dynamic system model including the bridle line system, the KCU and the tether. Because the focus of the present study is on the traction phase it can be assumed that the bridle lines and the tether are fully tensioned. Distributed loads such as gravity and aerodynamic drag will cause a certain deflection of line components. The degree of deflection depends on the length of the line and the ratio of the distributed load to the tensile force in the line. For the short bridle

5

lines deflection is only of minor importance, however, for long tethers the effect can be significant. According to [9], tether deflection can be neglected if the distance of the kite from the ground station does not exceed 100 m, which is the case for the present study. A schematic representation of the dynamic system model is shown in Fig. 5. The mass distribution of the wing is approximated by point

Table 1. Nominal parameters for the dynamic system model Parameter mA , mC m B , mD mKCU βAC βBD βAD , βBC dc ζc lt = rKCU LA , LC kc cc

zk

yk mA

are in fact very small. These numerical artifacts can be eliminated by imposing distance constraints between the point masses. The parameters values for the discussed model components are summarized in Table 1. The resulting dynamic system has 12 degrees of

mB

rk mC

xk mD

mKCU

Description discrete wing masses leading edge discrete wing masses trailing edge discrete mass KCU rel. damping coeff. A-C rel. damping coeff. B-D rel. damping coeff. A-D, B-C diameter bridle lines rel. damping coeff. bridle lines constant length tether rest lengths power lines stiffness bridle lines damping coefficient bridle lines

Fk

zw

freedom, assembled in the dynamic state vector y. The corresponding equations of motion are M y¨ = f, (12)

Fd m¨r

lt

m yw

1 2 FD,c

Fc

where M is the generalized mass matrix and f is the vector of the generalized forces. The equations of motion were symbolically derived in Matlab® .

Fg

xw

Figure 5. Dynamic system model with force components detailed for point mass mD attached to rear corner of left wing tip.

masses mA , mB , mC and mD positioned at the bridle line attachment points, taking into account that 70% of the total mass of the membrane structure is concentrated in the leading edge. The point mass mKCU represents the inertial properties of the control unit. The mechanical behavior of the bridle lines is modeled by linear spring-damper elements. The rest lengths LA and LC of the springs representing the power lines are constant, whereas the rest lengths LB and LD of the springs representing the steering lines can be modified individually to simulate steering, powering and de-powering of the wing. In the context of this paper, the steering lines are also labeled by their wing tip attachment as right and left steering lines, connecting to masses mB and mD , respectively. The forces acting on point mass mD are detailed in Fig. 5. The cable force Fc acting on mass mD is pointing towards mass mKCU and is computed as Fc = kc (lD − LD ) + cc l˙D , (10) where lD = krD − rKCU k is the actual length of the steering line, kc its stiffness and cc its damping coefficient. The force that the wing exerts on mass mD is calculated as Fk = −Fb with the reaction force given by the second block of equations in system (1). The aerodynamic drag force FD,c acts on the connected bridle line in the direction of the normal component va,n of the local apparent wind velocity. It is computed as FD,c =

Value 3.5 kg 1.5 kg 3 kg 0.02 0.002 0.015 0.0012 m 0.8 80 m 25 m 4.52 × 103 N/m 251.7 Ns/m

1 ρ CD lD dc v2a,n 2

(11)

where ρ denotes the air density, CD the aerodynamic drag coefficient of a cylinder in cross flow [36] and dc the diameter of the line. Half of the force is applied on mass mD and half of it on mass mKCU . Inertial effects have not been taken into account in the finite element model, which is why the gravitational force Fg and dynamic force m¨r acting on mass mD are explicitly included in the dynamic system model. The damping forces Fd are introduced between the wing tips to describe aerodynamic damping effects that are not covered by the aerodynamic load model. The damping constants are obtained by deriving the bending stiffness between the wing tips from numerical experiments and selecting appropriate relative damping coefficients βAC , βBD , βAD and βBC between the point masses to dampen the transverse wing tip relative motion. The stiff beams integrated in the wing tips connecting masses mA and mB as well as masses mC and mD can cause high-frequency oscillations in the dynamic model, although the corresponding displacements

D.

Numerical time integration

Equation (12) can be reformulated as y¨ = y¨ (y, y˙ , X, t) = M−1 f,

(13)

which indicates the functional dependence of the accelerations y¨ on the state vector y, the generalized velocities y˙ and the state vector X of the wing. The numerical integration over time is illustrated in Fig. 6. On input the FSI module receives the state vector X of the Initialization q0 , X0 Update flight state vector X X FSI problem

Flight controller y qb LB , LD

Fb Time integration y¨ (y, y˙ , X, t)

End time?

NO

New timestep ∆t

YES

Stop simulation

Figure 6. Numerical integration process of the iteratively coupled models.

6

wing and the displacements qb of the bridle attachment points. Once the finite element solver has converged to the static solution, it returns to the dynamic system model the reaction forces Fk = −Fb that act on the point masses mA , mB , mC and mD . Subsequently, the dynamic state vector y is advanced by a time step to the next discrete time of the system and the state vector X of the wing is updated. In order to fly the wing on a predefined trajectory and to keep the tension in the tether between a minimum and maximum value, the integration process includes steering and force controllers. Both implementations are described in the following section. The choice of the time integration scheme strongly affects the achievable accuracy, computational costs and stability of the simulation. Two different types of methods can be used to advance the state of the dynamic system model from one discrete time tn to the next tn+1 . Explicit methods calculate the state at tn+1 directly from state variables at tn , whereas implicit methods require solving a system of linear equations involving state variables at tn and tn+1 . As consequence, the computational cost per time step is low for explicit methods, however, the solution of stiff problems generally requires very small time steps because of the conditional numerical stability of the methods. Implicit methods, on the other hand, generally allow for large time steps because of their intrinsic stability. Another important aspect affecting the computational performance of the iterative solution process is the cost of calculating the static aeroelastic response of the wing per integration time step. A large time step generally implies larger changes in aerodynamic loading and shape deformation which has the consequence that the Newton-Raphson scheme employed by the FSI module will require more subiterations to advance the static equilibrium shape of the wing from tn to tn+1 . Numerical experiments have shown, that the large time steps of implicit integration methods lead to disproportionately expensive aeroelastic response calculations. Explicit methods, with smaller time steps and thus tighter coupling of model components, were able to achieve a higher overall performance. A fourth-order Runge-Kutta scheme, as implemented in the ODE45 Matlab® function, resulted in the shortest simulation times for the test cases considered in this work. Future developments might involve an adaptive implicit-explicit algorithm which can dynamically consider the time evolution of the driving aeroelastic forces.

III.

Flight control

The objective of the current study is to simulate the operation of a tethered wing on predefined flight trajectories typical for airborne wind energy systems. To achieve this, steering and force controllers - collectively denoted as flight controller - are included in the iterative coupling loop. The steering controller adjusts the actual length difference ∆l = lB −lD of the left and right steering lines such that the wing deforms asymmetrically and generates an aerodynamic turning moment. This controller emulates the functionality of the implemented kite control unit (KCU), as described in [7, 16, 17, 37]. The force controller adjusts the line lengths lB and lD equally, forcing the entire wing to pitch rotate (and deform symmetrically) in order to retain the generated aerodynamic force and resulting tether tension between certain minimum and maximum values. This force controller is not included in the implemented KCU but incorporated in the functionality of the winch controller of the ground station. It should be noted that the intention is not to develop an optimized control strategy for a kite power system, but to provide a combination of robust controllers to perform numerical simulations of actual flight maneuvers. A.

Steering controller

In order to analyze the steering of a kite it is suitable to decompose its flight motion into components along and perpendicular to the tether. When applying a spherical coordinate system, as illustrated in Fig. 7, this results in radial and tangential components of motion. The direction of the tangential component vk,τ is controlled by the KCU whereas the magnitude of the radial component vk,r is controlled by the ground station. In the current study the assumption of a constant tether length eliminates the radial component. Figure 7 also shows the control reference frame (xc , yc , zc ), which slightly differs from the kite-fixed reference. The zc -axis is in line with the tether, whereas the yc -axis is perpendicular, pointing towards the bridle attachment point on the leading edge of the

zc xc

zw

vw

zenith yc xs unit sphere

x∗s

ψ τ

ys δ

Rk

∗ R∗k ψ τ∗

h zs et

yw en h∗

R=1

O

xw

Figure 7. Mathematical concepts and definitions employed by the steering controller.

right wing tip. The control reference frame and the kite-fixed reference frame coincide only in the particular situation, where the bridle line attachment points A, B, C and D span a rectangle which is perpendicular to the tether. In order to non-dimensionalize the steering problem, the position rKCU of the steering unit is projected radially on the unit sphere. At this projected kite position Rk , the local tangential plane τ is spanned by axes xs and ys . The orientation of the kite relative to the wind reference frame is described by three Euler angles with an x − y0 − z00 rotation order. The third rotation is the yaw rotation of the kite and describes its heading h on the local tangential plane τ. The controller steers the kite such that the projected kite position Rk converges towards the desired trajectory on the unit sphere. Fig. 7 depicts this figure-eight trajectory, a Bernoulli lemniscate, which is a parametric curve that approximates well the power generating crosswind flight maneuvers [16, 17, 38]. The working principle of the controller is illustrated in Fig. 7 and can be summarized in five steps: (1) Find the point R∗k on the desired trajectory that has the shortest geodesic distance δ to point Rk . (2) At this point determine the vectors et and en that are tangential and normal to the desired trajectory and that span a local tangential plane τ∗ . (3) The new desired heading for the kite h∗ is represented as a linear combination of these vectors h∗ = βet + (1 − β)en .

(14)

The parameter β depends on the maximum permissible yaw rate of the kite, its current flight velocity and the distance δ to the trajectory, as explained in more detail in [16, 17]. (4) The heading angle ψ is measured from the xs -axis in the tangential plane τ and, similarly, the angle ψ∗ of the desired heading is measured from the x∗s -axis in the tangential plane τ∗ . (5) The error ε = ψ∗ − ψ is used as input for the controller. A proportional controller with gain K1 is introduced to control the desired length difference ∆l∗ between the steering lines, based on the error ε ∆l∗ = K1 (ψ∗ − ψ). (15) Two additional proportional controllers, both with gain K2 , are included to control the individual steering line lengths by controlling the rest lengths LB and LD of the right and left steering line springs based on the error ε0 = ∆l∗ − ∆l LB = LB,0 − K2 (∆l∗ − ∆l) LD = LD,0 + K2 (∆l∗ − ∆l)

(16) (17)

where LB,0 and LD,0 are the initial rest lengths of the springs. It should be mentioned that the described steering controller only resembles the more advanced tracking controller presented in [16, 17]. The ideas and the kinematic framework are the same, the terminology similar, however, the control laws differ substantially.

7

B.

Force controller

In order to avoid excessive tension forces in the steering lines, a simple proportional force controller with gain K3 is added to the steering controller, as depicted in Fig. 8. The objective of this controller is to keep the ratio between the forces in the steering and the power lines, µ=

FcB + FcD , FcA + FcC

(18) A.

at a desired value µ∗ = 0.3, which is characteristic for crosswind operation of the wing. The lengths of both steering lines are controlled simultaneously by adding the increment L3 = −K3 (µ∗ − µ)

(19)

to Eqns. (16) and (17). The interaction of the two controllers is illustrated in the block diagram of Fig. 8. + −

µ∗

0.34 which is due to the low gain employed in the force controller. The computations were performed on an Intel core i5 750 CPU with a clock speed of 2.8 GHz and a memory of 4 GB. Benchmark calculations with different numerical solvers showed that the explicit Runge-Kutta(4,5) algorithm with variable time step can achieve the shortest run times, about 27.5 times longer than the time duration of the simulated maneuver. The average integration time step was ∆t = 5.3 × 10−3 s.

Force controller

Steering mechanism

A central objective of the present study is the direct numerical simulation of the aeroelastic and flight dynamic response of a LEI tube kite to prescribed changes of the bridle geometry. According to [9], pulling on one of the steering lines not only increases the aerodynamic forces at the respective wing tip but also leads to a spanwise torsion of the entire wing. This asymmetric deformation significantly amplifies the aerodynamic yaw moment, forcing the wing to turn around the actuated tip. The mechanism is illustrated in Fig. 10 for pulling on the right steering line a

K3 Right tip ψ∗ + −

K1

∆l∗ + + −

1 K2

−1

+ −

LB,0 +

µ

LB

+ −

Kite +

LD

ψ

LD,0

Left tip A

FrL

B FrD

Yaw

Offset

FlL

C D

FlD

Steering controller

Figure 8. Block diagram of the steering and force controllers.

IV.

Results

The computational model and flight controller are implemented in Matlab® . The representative figure-eight flight maneuver selected for analysis is depicted in Fig. 9, showing the prescribed trajectory and the computed actual flight path. The kite is initially positioned in the

P0 P2 P1

P3 left turn

P4

right turn

Figure 9. Exterior view of the prescribed trajectory (dashed) and actual flight path (solid) with a total flight time of T = 20 s at a wind speed of vw = 5 m/s.

center of the figure-eight at point P0 and subsequently travels through points P1 , P2 , P3 and P4 to continue with the next figure. Accordingly, the kite performs a left turn during the first half of this maneuver and a right turn during the second half. The kite flies towards the ground in the side regions of the figure-eight, at |φ| ≈ φmax , and upwards in its center region, at φ ≈ 0. In practice, this downloop operation is used for equalizing the traction force during a figure-eight maneuver because the gravitational acceleration can partially compensate for the decreasing effective wind velocity in the side regions such that the kite can maintain a high flight speed [7]. The limitation of the basic proportional steering controller shows particularly in the lower altitude regions of the trajectory, from P1 to P2 and P3 to P4 , where the actual flight path deviates significantly from the prescribed trajectory. The force ratio µ varies between 0.26 and

Figure 10. Bottom view of a kite performing a right turn by pulling the right steering line (attached to point B).

The photo from below shows how the torsion of the wing leads to an increased angle of attack at the right wing tip, while decreasing the angle at the left tip. As consequence, the aerodynamic forces on the right tip increase, while decreasing at the left tip. However, the drag and the lift components contribute quite differently. The resulting yaw moment of the drag components, acting in the same direction, depends mainly on the force difference and is not significantly affected by the torsion of the wing. The lift components act in opposite directions and the resulting yaw moment increases linearly with the torsion-induced offset between the lines of action, as indicated in Fig. 10. The computational analysis presented in [9] shows that this torsion-based steering mechanism is responsible for the excellent maneuverability of the C-shaped kite. In the following, the computational simulation of a complete figureeight flight maneuver is interpreted on the basis of this mechanism. The deformation response of the wing model to a steering input is visualized in Fig. 11 applying magnified nodal displacements. The spanwise torsion is qualitatively similar to the behavior of the real wing depicted in Fig. 10. Also included is the distribution of nodal aerodynamic forces Fiaj , the resultant lift and drag components acting on the wing tips as well as the cable forces. Line forces and drag components are scaled by a factor of 2. Figure 12 shows the evolution of the local angle of attack for three representative wing sections. The left and right wing tips are characterized by angles αl and αr , respectively, whereas the midwing section is characterized by angle αm . The latter is directly related to the actual geometry of the bridle line system which is influenced by the force ratio µ defined in Eq. (18). The nearly constant value αm ≈ 7.5° proves the effectiveness of the force control strategy discussed in Sect. III. While flying the left turn during the first half of the figure-eight, the angle αl at the left tip is increased and the angle αr at the right tip is decreased. This situation inverses during the second half while the kite is flying a right turn. The increased angles of attack at the wing tips lead to increased aerodynamic drag and lift forces, as illustrated in Fig. 13. These results quantitatively support the steering mechanism proposed in [9] and discussed above. After passing point P1 and again after passing point P3 the lower one of the two tips angles temporarily drops to negative values and a In kite surfing references, left and right are generally defined from the perspective of the surfer on the ground, looking at the kite, whereas in this study the reversed definition from the kite perspective is adopted.

8

Aerodynamic force (N)

P0 600

left

P1

500 400

P2

right

P3

turn

P4

300 200 100 0 0

Yaw

turn

FLl FLr FDl FDr

10 Time (s)

5

15

20

FrL

Figure 13. Aerodynamic lift and drag force on the tips of the wing. FrD Right tip

FlD

P0 6

Offset

5

FcB FCc

FcA

Line force (kN)

FlL FcD Left tip

4 3

left

P1

turn

P2

right

P3

turn

P4

Ft FcA FcB FcC FcD

2 1

Figure 11. Perspective view of the twisted wing during a right turn.

Angle of attack (°)

P0 15

left

P1

turn

P2

right

P3

turn

0 0

5

10 Time (s)

15

20

P4

Figure 14. Tensile forces in tether and bridle lines.

10 5 0 -5 0

B. αl αr αm 5

10 Time (s)

15

20

Figure 12. Angle of attack of the left tip, right tip and middle section of the wing.

the corresponding lift component drops to low values. At the same time, both tip angles and the aerodynamic forces start to oscillate. This aeroelastic phenomenon, which can also be observed in flight tests, is caused by the severe wing torsion due to the strong increase of the steering input required to track the prescribed trajectory after passing points P1 and P3 (compare Fig. 9). As can be seen in Fig. 14, the resulting traction force of the wing varies between 6 kN, when passing the center of the figure-eight and 3 kN, when flying the downloop at the sides of the figure. It can be concluded, that the gravitational acceleration of the kite does not sufficiently compensate for the reduced effective wind speed at the sides of the figure-eight. The actual lengths lB and lD of the steering lines, the length difference ∆l = lB − lD and the steering input ∆L = LB − LD are displayed in Fig. 15 over the course of the flight maneuver. The variations of lB and lD are correlated and reflect the elastic response of the steering lines to the varying aerodynamic loading of the wing, superimposed by the adjustments of the rest lengths LB and LD by the force controller. Figure 16 shows the steering input ∆L and the flight dynamic response of the kite in terms of its yaw rate. In the upper half of the figure-eight, from P0 to P1 and P2 to P3 , the data is linearly correlated. In these parts of the maneuver, the kite is closely tracking the prescribed trajectory, while steering input and deformation of the wing are only moderate. A proportional relationship between yaw rate and steering input is also reported by [15–18], analyzing flight data from kites of different designs and ranging from ten to several hundred m2 surface area. The correlation deteriorates in the lower half of the figure-eight, from P1 to P2 and P3 to P4 , where the steering input is increased strongly to follow the prescribed trajectory. As discussed above, the resulting torsion of the wing is so pronounced that the angle of attack of one of the wing tips temporarily drops to negative values. This anomaly significantly affects the steering behavior of the wing.

Aerodynamic performance

The aerodynamic performance of the wing can be characterized by the integral lift to drag ratio L/D, which is plotted in Fig. 17. The ratio varies between a maximum value of L/D = 9, when the wing passes the center of the figure-eight at points P0 , P2 and P4 , and a minimum value of L/D = 7.5, when the wing flies from the sides of the figure-eight back to the center. This temporary drop in aerodynamic performance is mainly caused by the strong increase of drag at the wing tips as a result of the steering activities and the severe deformation. The high-frequency oscillations are related to the variations of the angle of attack which are generated by the force controller and which are visible in Fig. 12. As reported in [39], the measured lift to drag ratio of a similar type of kite is around L/D = 6, which is significantly below the range of values occurring in the simulation. Since the flight speed of a tethered wing in crosswind operation depends approximately linearly on L/D, the computed speeds of 18 m/s to 28 m/s are higher than actually achievable [5]. This can be explained as follows. Firstly, the aerodynamic model underestimates the effect of induced drag due to three-dimensional flow phenomena such as tip-vortices. Secondly, the geometric simplification of attaching the canopy to the center line of the leading edge tube instead of its top (see Fig. 4) artificially reduces the camber of the wing sections. Because this side effect decreases the aerodynamic drag more than the lift, the lift to drag ratio of the wing is increased. This should be corrected in a future revision of the structural model. Thirdly, the coarse discretization makes the finite element model stiffer than the real wing, which counteracts the deformation of the canopy and leads to a smaller camber than in reality. C.

Spanwise bending

The torsion of the wing in response to steering actuation is only one of the characteristic macro-scale deformation modes occurring during flight. Another one is the bending in spanwise direction resulting from the aerodynamic loading. This is illustrated in Fig. 18 by means of the transverse displacements of the bridle attachment points. The data indicates that the aerodynamic forces cause a large deformation of the tensile membrane structure from its initial unloaded shape. The heavily loaded center part of the wing is pulled outwards (in positive zk direction), tensioning bridle line system and tether and forcing the less loaded wing tips to move towards each other. This bending deformation, which essentially is an aerodynamic pretensioning of the structure, reduces the tip distance of the unloaded wing by an average value

9 P1

turn

P2

right

P3

turn

lB lD ∆l ∆L

25.5

P4 0.24 0.12

25.4

0.00

25.3

-0.12

25.2 0

10 Time (s)

5

P0 9.0 Lift to drag ratio (−)

left

Length difference (m)

Line length (m)

P0 25.6

-0.24 20

15

turn

P2

right

P3

turn

P4 1.5

0.2

1.0

0.2

0.5

0.0

0.0

-0.1

-0.5

-0.2 -0.3 0

∆L Yaw rate

-1.0 10 Time (s)

5

15

-1.5 20

Figure 16. Relation between steering input and yaw rate of the kite.

∆y = 1 m. This is about 20% of the design value given in Fig. 3, which clearly indicates the necessity of a geometrically non-linear modeling approach. Nearly symmetric tip deflections can be observed when the kite passes the center of the figure-eight, at points P0 , P2 and P4 , which corresponds to vanishing steering input ∆L and symmetric bridle geometry. Along the flight path, the transverse displacement of each wing tip varies in total by ∆y = 0.3 m, which is illustrated in Fig. 19 by two representative extreme shapes of the deformed wing. The variations of the tip displacements displayed in Fig. 18 do not show any notable oscillations. A similar behavior can be observed in flight tests. Even though an inflatable structure is very flexible, the aerodynamic pretensioning significantly increases its stiffness. During flight tests, wind gusts can inflict strong variations in wing loading, which can lead to large deformations, in extreme cases also to collapse of the tensile membrane structure. Gust loads are not included in the present model but would constitute an important extension to simulate the aeroelastic response of a tethered inflatable wing in a realistic wind field.

V.

Conclusions

In contrast to a rigid wing, a tethered inflatable wing is deformed substantially by the aerodynamic loading occurring during operation. Asymmetric changes in the bridle geometry directly affect the wing shape and can be employed efficiently for steering. The objective of the present study is to demonstrate the feasibility and performance of an approach for the coupled numerical simulation of the flight dynamics and structural deformation of a leading edge inflatable kite for traction power generation. Assuming a quasi-steady aeroelastic behavior, the wing is described by a static, geometrically nonlinear finite element model and combined with a steady aerodynamic load model. The reaction forces at the bridle line attachment points serve as input for a dynamic system model governing the flight dynamics of the wing. The results show that the computational model can reproduce the characteristic macro-scale bending and torsion deformation modes of a wing in a crosswind figure-eight flight maneuver. Compared to the unloaded wing, the deformations during flight are large, which justifies the use of a nonlinear structural model. Neglecting the inertial structural forces greatly reduces the computational effort, since the numerical time integration has to be performed only for the dynamic system model, which, compared to the finite element model of the wing has only a limited number of degrees of freedom. The mechanism of steering

P2

right

P3

turn

P4

7.5

10 Time (s)

5

15

20

Figure 17. Lift to drag ratio of the wing.

Transverse displacement (m)

P1

turn

8.0

7.0 0

Yaw rate (rad/s)

Steering input (m)

left

P1

8.5

Figure 15. Actual lengths of steering lines, lenght difference and steering input. P0 0.3

left

P0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 0

left

P1

turn

P2

right

P3

turn

P4

C D

B A 5

10 Time (s)

15

20

Figure 18. Transverse (yk -direction) displacements of the bridle line attachment points.

is analyzed quantitatively and the results confirm the importance of wing torsion. For moderate degrees of deformation a linear relation between steering input and yaw rate is identified, which is consistent with previous experimental studies. The limiting factor in the computational approach is the accuracy of the aerodynamic wing loading model. The present implementation does not account for three-dimensional flow effects and employs an approximate method to redistribute the aerodynamic force on a wing section along its surface. It can thus be expected that the prediction quality can be increased substantially by addressing the deficiencies in the wing loading model.

VI.

Acknowledgments

This research was supported financially by the Rotterdam Climate Initiative. The authors would like to thank Wubbo J. Ockels for his pioneering contributions to the airborne wind energy research at Delft University of Technology and for inspiring and supporting the various development activities. Much appreciated are the contributions of former colleagues Jeroen Breukels, for valuable feedback, Claudius Jehle, for support during the controller development and Arend L. Schwab for creative input during the conceptual phase.

References [1] Archer, C. and Caldeira, K., “Global assessment of high-altitude wind power,” Energies, Vol. 2, No. 2, 2009, pp. 307–319. doi:10.3390/en20200307. [2] Jacobson, M. Z. and Archer, C. L., “Saturation wind power potential and its implications for wind energy,” Proceedings of the National Academy of Sciences, Sept. 2012. doi:10.1073/pnas.1208993109. [3] Archer, C. L., Monache, L. D., and Rife, D. L., “Airborne wind energy: Optimal locations and variability,” Renewable Energy, Vol. 64, 2014, pp. 180– 186. doi:http://dx.doi.org/10.1016/j.renene.2013.10.044. [4] Fagiano, L. and Milanese, M., “Airborne Wind Energy : an overview,” Proceedings of the 2012 American Control Conference, IEEE, Montr´eal, Canada, June 2012, pp. 3132–3143. [5] Schmehl, R., Noom, M., and van der Vlugt, R., “Traction Power Generation with Tethered Wings,” Airborne Wind Energy, edited by U. Ahrens, M. Diehl, and R. Schmehl, Springer, Berlin-Heidelberg, 2013, pp. 23–45. doi:10.1007/978-3-642-39965-7 2.

10 [19] Baayen, J. H., “Vortexje - an open-source panel method for co-simulation,” 2012, http://arxiv.org/abs/1210.6956v3. [20] Gohl, F. and Luchsinger, R. H., “Simulation Based Wing Design for Kite Power,” Airborne Wind Energy, edited by U. Ahrens, M. Diehl, and R. Schmehl, Springer, Berlin-Heidelberg, 2013, pp. 325–338. doi:10.1007/978-3-642-39965-7 18. zk

yk

Figure 19. Frontview of the spanwise bending of the wing.

[6] Luchsinger, R. H., “Pumping Cycle Kite Power,” Airborne Wind Energy, edited by U. Ahrens, M. Diehl, and R. Schmehl, Springer, BerlinHeidelberg, 2013, pp. 47–64. doi:10.1007/978-3-642-39965-7 3. [7] van der Vlugt, R., Peschel, J., and Schmehl, R., “Design and Experimental Characterization of a Pumping Kite Power System,” Airborne Wind Energy, edited by U. Ahrens, M. Diehl, and R. Schmehl, Springer, BerlinHeidelberg, 2013, pp. 403–425. doi:10.1007/978-3-642-39965-7 23. [8] Fechner, U. and Schmehl, R., “Design of a Distributed Kite Power Control System,” Proceedings of the 2012 IEEE International Conference on Control Applications, IEEE, Dubrovnik, Croatia, 2012, pp. 800–805. doi:10.1109/CCA.2012.6402695. [9] Breukels, J., An Engineering Methodology for Kite Design, Ph.D. thesis, Delft University of Technology, 2011, http://resolver.tudelft.nl/uuid: cdece38a-1f13-47cc-b277-ed64fdda7cdf. [10] J.P. Meijaard, W.J. Ockels, A. S., “Modelling of the dynamic behaviour of a laddermill, a novel concept to exploit wind energy,” 3rd International Symposium on Cable Dynamics, Trondheim, Norway, August 1999, pp. 229–234. [11] Williams, P., Lansdorp, B., and Ockels, W. J., “Modeling, Simulation, and Testing of Surf Kites for Power Generation,” AIAA-Paper 2008-6693, 2008. doi:10.2514/6.2008-6693. [12] Terink, E. J., Breukels, J., Schmehl, R., and Ockels, W. J., “Flight Dynamics and Stability of a Tethered Inflatable Kiteplane,” AIAA Journal of Aircraft, Vol. 48, No. 2, 2011, pp. 503–513. doi:10.2514/1.C031108. [13] Dadd, G. M., Hudson, D. A., and Shenoi, R. A., “Determination of kite forces using three-dimensional flight trajectories for ship propulsion,” Renewable Energy, Vol. 36, 2011, pp. 2667–2678. doi:10.1016/j.renene.2011.01.027. [14] de Groot, S. G., Breukels, J., Schmehl, R., and Ockels, W. J., “Modeling Kite Flight Dynamics Using a Multibody Reduction Approach,” AIAA Journal of Guidance, Control and Dynamics, Vol. 34, No. 6, 2011, pp. 1671– 1682. doi:10.2514/1.52686. [15] Erhard, M. and Strauch, H., “Control of Towing Kites for Seagoing Vessels,” IEEE Transactions on Control Systems Technology, 2012. doi:10.1109/TCST.2012.2221093. [16] Jehle, C., Automatic Flight Control of Tethered Kites for Power Generation, Master’s thesis, TU Munich, 2012, http://mediatum.ub.tum.de/?id= 1185997. [17] Jehle, C. and Schmehl, R., “Applied Tracking Control for Kite Power Systems,” Accepted for publication in AIAA Journal of Guidance, Control and Dynamics, 2013. [18] Fagiano, L., Zgraggen, A., Morari, M., and Khammash, M., “Automatic Crosswind Flight of Tethered Wings for Airborne Wind Energy: Modeling, Control Design, and Experimental Results,” IEEE Transactions on Control Systems Technology, Vol. PP, No. 99, 2013, pp. 1–1. doi:10.1109/TCST.2013.2279592.

[21] Furey, A. and Harvey, I., “Evolution of neural networks for active control of tethered airfoils,” Advances in Artificial Life, 2007, pp. 746–755. doi:10.1007/978-3-540-74913-4 75. [22] Williams, P., Lansdorp, B., and Ockels, W. J., “Flexible Tethered Kite with Moveable Attachment Points, Part I: Dynamics and Control,” AIAA-Paper 2007-6628, 2007. doi:10.2514/6.2007-6628. [23] Breukels, J., Schmehl, R., and Ockels, W., “Aeroelastic Simulation of Flexible Membrane Wings based on Multibody System Dynamics,” Airborne Wind Energy, edited by U. Ahrens, M. Diehl, and R. Schmehl, Springer, Berlin-Heidelberg, 2013, pp. 287–305. doi:10.1007/978-3-642-399657 16. [24] Schwoll, J., Finite Element approach for statically loaded inflatable kite structures, Master’s thesis, Delft University of Technology, 2012. [25] Chatzikonstantinou, T., “Numerical analysis of three-dimensional non rigid wings,” AIAA-Paper 89-0907, 1989. doi:10.2514/6.1989-907. [26] M¨uller, S., Modellierung, Stabilit¨at und Dynamik von Gleitschirmsystemen, Phd thesis, TU Munich, 2002. [27] Bosch, H. A., Finite Element Analysis of a Kite for Power Generation, Master’s thesis, Delft University of Technology, 2012, http://resolver.tudelft.nl/ uuid:888fe64a-b101-438c-aa6f-8a0b34603f8e. [28] Bosch, A., Schmehl, R., Tiso, P., and Rixen, D., “Nonlinear Aeroelasticity, Flight Dynamics and Control of a Flexible Membrane Traction Kite,” Airborne Wind Energy, edited by U. Ahrens, M. Diehl, and R. Schmehl, Springer, Berlin-Heidelberg, 2013, pp. 307–323. doi:10.1007/978-3-64239965-7 17. [29] Breukels, J. and Ockels, W. J., “A Multi-Body System Approach to the Simulation of Flexible Membrane Airfoils,” Aerotecnica Missili Spazio, Vol. 89, No. 3, 2010, pp. 119–134. [30] Allman, D. J., “Evaluation of the constant strain triangle with drilling rotations,” International Journal for Numerical Methods in Engineering, Vol. 26, No. 12, 1988, pp. 2645–2655. doi:10.1002/nme.1620261205. [31] Felippa, C., “A study of optimal membrane triangles with drilling freedoms,” Tech. Rep. CU-CAS-03-02, Center for Aerospace Structures, College of Engineering, University of Colorado, 2003. [32] Allman, D. J., “A simple cubic displacement element for plate bending,” International Journal for Numerical Methods in Engineering, Vol. 10, No. 2, 1976, pp. 263–281. doi:10.1002/nme.1620100204. [33] Reddy, J. N., An Introduction to Nonlinear Finite Element Analysis, Oxford University Press, USA, 2004. [34] Crisfield, M. A., Non-Linear Finite Element Analysis of Solids and Structures, Wiley, 1996. [35] Spierenburg, G., Continued development and experimental validation of a kite design and simulation tool, Master’s thesis, Delft University of Technology, 2005. [36] Hoerner, S., “Fluid dynamic drag,” USA: Hoerner Fluid Dynamics, 1965. [37] Baayen, J. H. and Ockels, W., “Tracking control with adaption of kites,” IET Control Theory & Applications, Vol. 6, No. 2, 2012, pp. 182–191. doi:10.1049/iet-cta.2011.0037. [38] Ruiterkamp, R. and Sieberling, S., “Description and Preliminary Test Results of a Six Degrees of Freedom RigidWing Pumping System,” Airborne Wind Energy, edited by U. Ahrens, M. Diehl, and R. Schmehl, Springer, Berlin-Heidelberg, 2013, pp. 443–458. doi:10.1007/978-3-642-399657 26. [39] Ruppert, M., Development and validation of a real time pumping kite model, Master’s thesis, TU Delft, 2012.