Kite as a beam: A fast method to get the flying shape

Designing new large kite wings requires engineering tools that can account for flow-structure interaction. Although a fully coupled simulation of deformable membrane structures under aerodynamic load is already possible using Finite Element and Computational Fluid Dynamics methods this approach is computationally demanding. The core idea of the present study is to approximate a leading edge inflatable tube kite by an assembly of equivalent beam elements. In spanwise direction the wing is partitioned into several elementary cells, each consisting of a leading edge segment, two lateral inflatable battens, and the corresponding portion of canopy. The mechanical properties of an elementary cell—axial, transverse shear, bending, and torsion stiffness—and the chordwise centroid position are determined from the response to several imposed elementary displacements at its boundary, in the case of a cell under an uniform pressure loading. For this purpose the cells are supported at their four corners and different non-linear finite element analyses and linear perturbation computations are carried out. The complete kite is represented as an assembly of equivalent beams connected with rigid bodies. Coupled with a 3D non-linear lifting line method to determine the aerodynamics this structural model should allow predicting the flying shape and performance of new wing designs. Alain de Solminihac · Alain Nême ( ) · Chloé Duport · Jean-Baptiste Leroux · Kostia Roncin · Christian Jochum ENSTA Bretagne, FRE CNRS 3744, IRDL, 29200 Brest, France e-mail: alain.neme@ensta-bretagne.fr Yves Parlier Beyond the sea, 1010 avenue de l’Europe, 33260 La Teste de Buch, France e-mail: yves.parlier@beyond-the-sea.com


Introduction
The aim of this research project is the development of tethered kite systems as auxiliary devices for the propulsion of merchant ships. It is part of the beyond the sea R program which is lead by the Institut de Recherche Dupuy de Lôme (IRDL) of EN-STA Bretagne. The goal is to design leading edge inflatable tube kites with surface area larger than 300 m 2 . This requires a significant upscaling of common sports kites which generally do not exceed a surface area of 30 m 2 . This upscaling process raises several issues: What are the relevant physical effects to take into account? Is it possible to use the same materials as for sports kites? Which geometry should the bridle system have?
Point mass and rigid body models have been used for real-time or faster-thanreal-time simulation of the kite dynamics [7,9]. A typical application of these type of models is control engineering or flight path optimization. However, to get a deeper understanding of the steering behavior and aerodynamic performance of a highly flexible wing the shape deformation plays a crucial role. Breukels [3,4] developed an engineering model of a deformable flying kite, discretizing the tubular frame by chains of rigid bodies connected by rotational springs and the canopy by arrays of elastic springs and damper elements. All mechanical properties were derived from basic experiments and the aerodynamic load distribution was prescribed by an empirically determined correlation framework. The approach allows modeling of aeroelastic effects.
Bosch [2] applied a geometrically non-linear finite element framework to the kite, discretizing the tubular frame by beam elements and the canopy by custom-made shell elements. This model was used to determine the quasi-static deformation resulting from changes in the boundary conditions, such as aerodynamic loading and steering line displacements. However, only macro-scale Fluid-Structure Interaction (FSI) effects of spanwise torsion and bending of the wings were taken into account. Gaunaa [8] developed a computationally efficient method for determining the aerodynamic performance of kites. The approach iteratively couples a Vortex Lattice Method (VLM) with 2D airfoil data to account for the effects of airfoil thickness and of viscosity. Deformation of the wing is not considered.
The aim of the present study is to develop an engineering tool which enables kite designers to efficiently determine the flying shape of new kites. Given a few design parameters such as the global wing shape, the material used and the wind conditions, it should be possible to predict the flying shape and aerodynamic performance. Structural non-linearity and macro-scale FSI calculations are conducted as major influence factors. Contrary to the point mass and rigid body approaches, the proposed method underlines the importance of considering the inflatable kite as a deformable membrane structure. The method can be used to identify critical aerodynamic peak loads and to take design measures to alleviate these.
In Sect. 4.2 and 4.3 the wing design and basic methods are first outlined, introducing the dicretization concept of the elementary cell, then followed by the identification of mechanical and inertial properties and the assembly of several elementary cells into a model of the complete wing. In Sect. 4.4 and 4.5 results are presented, discussed and interpreted, elaborating also on the fast potential flow-based method used to derive the instantaneous aerodynamic loading of the wing. The preliminary content of the present chapter has been presented at the Airborne Wind Energy Conference 2015 [17].

General Design Parameters
In this section the general problem definition is outlined. The kite design is developed in several steps starting from a 3D baseline. The required material properties are then discussed, followed by a specification how the aerodynamic loading is determined for different operational modes of the kite.

Design Geometry
The spanwise shape of the wing design is defined by a 3D curve. The chord, sweep, dihedral, and twist of the wing are specified by evolution laws along this baseline. The inflatable tubular frame is detailed by specifying the attachment points of the inflatable battens at the leading edge tube as well as all tube cross section geometries. Each pair of neighboring battens and the corresponding part of the leading edge tube spans a wing section. To complete the definition of the design geometry the design camber of each of these wing sections is defined. This property describes the maximum deviation of the canopy from the mean chord of the wing section. An example of the resulting wireframe representation of the wing is shown in Fig. 4.1.

Baseline
Trailing edge

Material Properties
The fabric mechanical properties for the inflatable beams and the canopy are defined using some data per unit mass (J/kg) for specifying the specific Young's modulus E m , and some data per unit area (kg/m 2 ) for specifying the fabric density µ. The Poisson's ratio ν of the fabric should also be specified. In this study only isotropic materials are considered, but the present model could be extended in the case of anisotropic materials.

Relative Wind Conditions
The relative flow conditions at the wing are defined by the apparent wind velocity with v w denoting the true wind velocity and v p the ship velocity, both known properties, and v k denoting the kite velocity relative to the ship. Kites can be used in two different flight modes to generate a traction force for towing ships. In static flight mode the kite has a fixed position with respect to the ship and the apparent wind velocity can be readily calculated from Eq. (4.1) by setting v k = 0. In dynamic flight mode the kite is operated perpendicularly to the tether and the kite velocity is a variable. It is possible to use a simple dynamic flight model, such as the zero mass model [6,[12][13][14], to calculate v k as a function of time and to use this in Eq. (4.1) to derive the apparent wind velocity.
The model described further in the next section requires an a priori estimation of the pressure loading of the canopy, because its geometrical stiffness must be considered. From the apparent wind velocity, and given an estimate a priori of the aerodynamic lift coefficient, this can be achieved by calculating with ρ denoting the air density and C L the aerodynamic lift coefficient. Given the relatively high lift-to-drag ratio of the wings involved, and given an approximate pressure loading is only required, the effect of the aerodynamic drag coefficient is neglected here.

Structural Model of the Wing
In this section the structural model of the wing is built up in steps, starting from individual elementary cells which are assembled into a structural model of the entire flexible wing.

Elementary Cell Concept
The structural discretization of the kite is based on a spanwise division of the wing into sections. The proposed concept of the elementary cell accounts for the particular structural design of a membrane wing with inflatable tubular frame. As illustrated in Fig. 4.2 each cell is composed of a segment of the inflatable leading edge, two inflatable battens and the corresponding portion of the canopy. The mechanical be- Because the geometry of the wing is double-curved an elementary cell Q ′ : To further simplify the cell geometry the planar approximation Q : L L , L R , T R , T L is introduced. Using the midpoints M 1 and M 2 on the left and right batten segments L ′ L , T ′ L and L ′ R , T ′ R the spanwise dimension L and the mean chord H of the planar approximation Q are defined as A local coordinate system (e 1 , e 2 , e 3 ) is defined by the unit vector along the spanwise direction the unit vector perpendicular to the plane and a third unit vector defined as cross product

Equivalent Beam Concept
The equivalent beam is introduced to describe the mechanical behavior of the elementary cell by means of an idealized structural object. The following beam properties are identified on the basis of finite element analysis of the elementary cell under various loads: • Beam centroid distance from the leading edge, • Tension/Compression stiffness, • Bending stiffness, • Torsion stiffness, • Shear coefficients.
The structural analysis is performed with the finite element solver Abaqus TM .

Finite Element Model of the Elementary Cell
As a conclusion of a convergence analysis the canopy of the elementary cell is discretized by 2000 rectangular linear membrane elements. The mechanical properties used for the canopy are the in-plane stiffness E C = µ C E m,C and the Poisson ratio ν C . The subscript C indicates properties of the canopy. It is possible to adapt these mechanical properties for the different regions of the canopy, as for instance at the trailing edge if canopy reinforcement effects have to be investigated. The canopy of the elementary cell is supported by the leading edge tube and two battens. These inflatable elements are modeled as straight beams and discretized by 200 linear beam elements in total for three tubes. Starting from the known beam radius R, fabric stiffness E B = µ B E m,B and Poisson ratio ν B , where subscript B indicates properties of the beam, the section properties are estimated as: • Elongation stiffness: 2πRE B , • Bending stiffness: πR 3 E B , • Transverse shear stiffness [5] : • Torsion stiffness: Because in the final wing model the elementary cells are connected the stiffness of the finite element beam representing a batten is only 50% of the stiffness of the full batten. Underlying is a linear superposition assumption. In this study, it is assumed that the kite design geometry allows to consider the tips of the wing as battens. For these, the stiffness of the corresponding finite element beam is 100% of the stiffness of the full batten.

Pressurization of the Elementary Cell
The geometrical stiffness of the canopy must be considered because it is comparable to the stiffness of the beam frame. The initial shape of the canopy before applying the pressure loading is expressed in the Cartesian frame (e 1 , e 2 , e 3 ) with origin at L L with x 3 given by the following analytic expression and λ denoting the design camber of the canopy with a value of λ ≈ 5%. The first computation step is a non-linear geometrical analysis. The four corners (T L , L L , L R , T R ) are clamped in space and the elementary cell is loaded with the estimated homogeneous pressure as described in Sect. 4.2.3.
Since membrane elements have no bending stiffness, a damping factor of 5 × 10 6 is introduced in the Abaqus TM simulation [16] to achieve convergence of the nodal force balance at the end of the time step (100 seconds). Then a second computation step is conducted without damping to check the validity of the obtained solution. A representative simulation result is shown in Fig. 4.3. As a last step the characteristics of the elementary cell under homogeneous pressure loading are determined.

Computation in Linear Perturbation Mode
Starting from this pressurized structure, five linear perturbation calculation cases are completed in order to evaluate the stiffnesses of the elementary cell with respect to the different global degrees of freedom. The cases are listed in Table 4  Numerical results depend linearly on a since a linear perturbation mode is used. Reaction forces at the right corner points, T R and L R , are measured for each load case in the direction of the elementary displacement. F TR,X is the reaction force at the trailing edge and F LR,X is the reaction force at the leading edge for load case X.
The computed deformation of the elementary cell is shown for two representative load cases.

Identification of Beam Properties
The objective of the equivalent beam is to model the mechanical behavior of the elementary cell. The specific load cases (a), (b), (c), (d) and (e) are studied. Static equilibrium on a Timoshenko beam, without warping effect, allows getting relationships between measured reaction forces, elementary displacements and equivalent beam properties. Solving this system of equations, the equivalent beam can be completely described. For better understanding, this approach is presented for the torsion load case and the equivalent beam extremities B 1 and B 2 are given by The stretching stiffness is calculated as The in-plane bending stiffness about e 3 is determined as The strain energy ratio of transverse shear stiffness along e 3 and bending stiffness along e 2 is conventionally evaluated as 12EI 2 /(GA 03 L 2 ). For all the studied cases, this ratio is approximately 3, which is expected for standard leading edge inflatable tube kites. According to this property the transverse shear stiffness along e 3 can be estimated as (4.21) The transverse shear stiffness along e 2 is given by and the bending stiffness about e 2 is evaluated as

Wing Assembly
To build the kite structure, equivalent beams representing elementary cells are gathered and connected together with rigid bodies. This method is illustrated in Fig. 4.6. The equivalent beams edges are not at the same position for two successive elementary cells. It is, then, assumed here that two neighbor beams share similar displacements and rotations at their extremities similarly to a virtual rigid body connecting these two sections.

Case of Study
The data used for the case of study is summarized in

Flow Model-Non-Linear 3D Lifting Line Method
The non-linear 3D lifting line is based on an extension of Prandtl's lifting line theory. This extension is intended for wings with variable dihedral and sweep angles. Leloup introduces in [14] a linear implementation while the present method is considering the non-linearity of the aerodynamic lift coefficient. The finite wing and its wake are represented by a set of horseshoe vortices of different circulation strengths Γ . The aim of the algorithm presented below is to calculate the circulation of each horseshoe vortex. Once obtained, the local effective flow for each wing section allows calculating the local aerodynamic forces and torques along the wing span. The numerical iterative solution is based on Anderson [1, Chap. 5, Sect. 5.4], the calculation of effective local incidence angles was adapted to the cases of wings which are non-straight and non-planar. The horseshoe vortices used for discretisation and calculation of their influences are derived from Katz and Plotkin [10, Chap. 12, Fig.  12.2 (a)]. The wing is divided in a finite number of parallel sections, each one represented by a horseshoe vortex. A horseshoe vortex consists of six vortex segments. The bound vortex is located at the quarter chord length, carefully perpendicular to the plane of the considered section. Each of the two trailing vortices are separated into two parts: the first one extends parallel to the chord over one chord length and the second one extends parallel to the local free stream over several chord lengths. Finally the starting vortex closes the horseshoe. It is important to note that even with a swept wing, the bound vortex along the lifting line is orthogonal to the two adjacent trailing vortices. This is illustrated in Fig. 4.7. This leads to a piecewise constant discretisation of the lifting line, but it is necessary to have a correct match between the local lift calculated from the Kutta formula or from the polar of the section.

Horseshoe vortices
Leading edge Trailing edge Local aerodynamic forces x z y v ∞ α In order to calculate the aerodynamic forces, the circulation is initialized by an elliptic distribution along the wing span. With the Biot-Savart law, the induced velocities by each vortex segment can be calculated and summed at each point of the lifting line. Combined with the local free stream velocity, the effective wind and the effective incidence angle are obtained for each section. The current local bound vortex strength is then calculated using the polar of the section, which leads to the local lift force, and via the Kutta formula, which converts this force in local circulation. The circulation value is ultimately updated by weighting between current and previous values using a damping factor. This whole process is repeated until the circulation distribution converges. The lift, drag and torque of each section of the wing are then post processed with the converged local circulation value, which leads to integrated local loads, and after being carefully summed, to global loads of the wing.

Abaqus TM Procedure
The structural analysis is conducted using the commercial solver Abaqus TM . The computation of the equivalent beam deformation under aerodynamic loading is performed with a large-displacement formulation from the initial configuration of the equivalent beam which accounts for its stress-free geometry. The large-displacement formulation of Timoshenko beam elements used in Abaqus TM [16] is based on a multiplicative decomposition of the deformation gradient into a stretch part (F s ) and a distorsion part (F d ). The strain tensor is obtained by addition of the logarithm of F s and the Green-Lagrange formula applied to F d . No artificial damping forces were introduced into the finite element model. Since the geometrical location of the finite element beam lies on the lifting line, its local section direction n 1 is determined with the orthogonal projection of the point M, defined as the geometric center of the beam element, on the equivalent beam which is located at the distance D (see Eq. (4.15)) from the leading edge. If P represents the projection of M, it can be determined from If t stands for the unit vector along the beam element axis, the unit vector n 1 is obtained from (4.25) The second local section direction of the beam element n 2 is such that We assume that the location of the beam element section centroid is expressed in the local beam element frame (t, n 1 , n 2 ) as The beam element section properties are the same as in the Eqs. (4.18) to (4.23) assuming the local beam element frame (t, n 1 , n 2 ) is matching the equivalent beam frame (e 1 , e 2 , e 3 ).

Boundary and Wind Conditions
To model the closing and opening of the kite under load the specific boundary conditions listed in Table 4.3 are chosen. By definition the apparent wind velocity is aligned with the x-axis as illustrated in Fig. 4.1. It has a value of 30 m/s at an air density of 1.2kg/m 3 . No twist is considered for the stress-free geometry of the kite and the wind is parallel to its symmetry

Fluid-Structure Coupling
An iterative algorithm [15] is used in a single artificial time increment corresponding to the kite evolution from its stress-free configuration up to its deformed configuration under aerodynamic loading. The first beam loading is computed with the 3D lifting line method considering the stress-free configuration. A similar procedure is used in [18]. The same line is used for both flow model and structure calculation and both lines have the same mesh. Fluid computations provide nodal aerodynamic forces and moments reduction whereas solid calculations determine nodal displacements and rotations. Note that the deformation of the kite does not change the 2D characteristics of the wing section used for fluid computations. The convergence of the procedure is observed through two physical values: lift and kite closing.

Results
The CPU times observed on a classical computer 1 are respectively 0.12 s and 1.3 s to obtain the non-linear lifting line and Abaqus TM converged solutions. Generally, six fluid-structure coupling loops are required to achieve the convergence, as can be seen in Fig. 4.8. This is quite similar to the convergence observed on former study with a shell finite element modeling of the canopy [11]. In Fig. 4.9, the torsion of the canopy can be observed whereas in Fig. 4.10 the opening is shown. The presented design does not contain any bridle system hence such large displacements can be noticed. +0.000e+00 +3.669e-01 +7.339e-01 +1.101e+00 +1.468e+00 +1.835e+00 +2.202e+00 +2.569e+00 +2.936e+00 +3.302e+00 +3.669e+00 +4.036e+00 +4.403e+00 Deformed state Initial state

Discussion
This method allows a first estimation of the flying shape, the drag and the lift of a deformable kite. As highlighted in Fig. 4.8 the global lift force is around 40% lower for a soft kite (converged point) than for a rigid kite (first point). This result tends to justify the soft kite approach for the simulation of kite performances. For realistic displacements, the bridle system and tethers must also be taken into account. It is +0.000e+00 +3.669e-01 +7.339e-01 +1.101e+00 +1.468e+00 +1.835e+00 +2.202e+00 +2.569e+00 +2.936e+00 +3.302e+00 +3.669e+00 +4.036e+00 Deformed state Initial state

Conclusion
In this study the complex structural behavior of a soft kite was simplified to a simple arrangement of beams. The parameters of the beams were calculated from finite element analysis of so-called elementary cells, which model the canopy of a single kite cell, under homogeneous pressure. This pressure was derived from the global lift coefficient of the initial kite geometry.
Coupled with a fluid model, the simplified structure model approach presented in this study allows a prediction of the flying shape and helps obtaining a better understanding of the main phenomena which have to be considered. It is as well a quick (a couple of minutes) and convenient way to get a first estimation of the kite performance accounting for fluid-structure interaction.
However, the "kite as a beam" model has not been compared to more detailed structural models. This analysis is currently ongoing. Additionally, the "kite as a beam" approach presented here does not directly address local aspects like stresses and strains in the canopy and in the inflatable structure. These aspects have to be investigated with fully coupled FEA / CFD computations. Overall validation requires relevant experiments that are currently under progress at the institute.
The next step to extend the "kite as a beam" model would be the inclusion of bridles and tethers to improve their design and for better towing force estimations. In parallel, it will be necessary to develop a more realistic beam frame model for the kite structure. A parametric formulation of the influence of the geometry of the canopy on the kite stiffness will be developed. Hence, the stiffness of the elementary cells will depend on the aerodynamic pressure. The new model should also enable an improvement of the design of the inflatable leading edge and the battens according to stress limit and buckling condition.