leurs bonnes

Abstract

The present thesis is part of the beyond the sea® project which aims to develop tethered kite systems as auxiliary devices for ship propulsion in order to reduce greenhouse gas emissions and fuel costs. In a design stage, it is necessary to be able to estimate the traction that can provide a kite according to its trajectory. As a kite is a flexible structure, its geometry depends on its aerodynamic loading, which depends on its geometry. For that reason, fluid-structure interaction has to be taken into account to calculate the flying shape and aerodynamic performances of the wing. Knowing the computational time needed to carry out a fully coupled simulation using Finite Element and Computational Fluid Dynamics methods; it can be very useful in a design stage to use fast, reliable and simple models to estimate rapidly the kite performances. Moreover fast simulations allow the studies of different trajectories and the determination of the critical loading case, where more complex and more time-consuming models can be used.

A 3D Non-Linear Lifting Line model has been developed, based on the lifting line of Prandtl. This extension of the Prandtl lifting line is intended to deal with nonstraight kite wings, with dihedral and sweep angles variable along the span and take into account the non-linearity of the lift coefficient. The model is also able to deal with any general flight kinematic including translation velocities and turning rates, diverse angles of incidence or angles of sideslip.

The model has been checked with 3D RANSE simulations over various geometries and produces satisfactory results for range of incidence and sideslip up to 15°, with typical relative differences of few percent for the overall lift. The local results, needed for further fluid-structure interaction, are also correctly estimated; the model is able to predict the position of the minimum and maximum loading along the span, even in sideslip. Special attention has been given to the estimation of the accuracy of the provided numerical results.

The model has also been introduced in a new equilibrium iterative procedure for the kite, which has been used to estimate the loading of a dynamic kite flying on a circular path. This case points out the influence of the rotation rates on the loading of the kite. A parametric study has finally been conducted to analyze the evolution of the aerodynamic coefficients and the flight characteristics according to the tether length or the radius of the trajectory.

Simultaneously, a structure model has been developed. The core idea of the Kite as a Beam model is to approximate a Leading Edge Inflatable kite by an assembly of equivalent beam elements. Over the span, the kite is divided into several elementary cells, each one composed of a portion of the inflatable leading edge, two inflatable battens and the corresponding canopy. First of all, the cell is put under pressure and then subject to different linear displacement perturbations in order to calculate the tangent stiffness of the equivalent beam.

iii The Kite as a Beam model has been compared to a complete kite Finite Element model over elementary comparison cases. The structure models are put under a constant pressure and then the same elementary displacement is applied, a translation along one of the three global axes or a rotation of the wing tip. The results shows the behavior differences of the two models, for example the torsion stiffness is globally overestimated by the Kite as a Beam model.

Eventually, the Kite as a Beam model coupled with the 3D Non-Linear Lifting Line model is compared to the complete Finite Element model coupled with the 3D Non-Linear Lifting Line model. The gain in computation time is really significant as the Kite as a Beam model takes only a few minutes to compute while the complete Finite Element model needs one or two hours. However, the results show the necessity of model calibration to make both models coincide.

Résumé

La présente thèse s'inscrit dans le cadre du projet beyond the sea®, projet qui a pour but de développer l'exploitation de cerfs-volants pour la propulsion auxiliaire de navires. En effet, afin de réduire l'émission de gaz à effet de serre, la consommation d'énergie fossile (et donc le coût du transport maritime), l'utilisation de cerfs-volants pour la traction des navires apparaît comme une solution réalisable [START_REF] Wellicome | Ship propulsive kites: an initial study[END_REF]; [START_REF] Naaijen | On the power savings by an auxiliary kite propulsion system[END_REF]; [START_REF] Leloup | A continuous and analytical modeling for kites as auxiliary propulsion devoted to merchant ships, including fuel saving estimation[END_REF]; [START_REF] Erhard | Control of Towing Kite for Seagoing Vessels[END_REF]). Le projet beyond the sea®, fondé en 2007, a reçu le soutien de l'Agence de l'Environnement et de la Maîtrise de l'Énergie en 2014 et associe des entreprises privées et des partenaires académiques. L'utilisation de cerfs-volants présente de nombreux avantages, le système d'accroche et de contrôle du kite prend moins de place sur le pont d'un navire qu'un mât et ses voiles par exemple, un cerf-volant peut aller chercher le vent en altitude, où il est plus régulier et plus fort [START_REF] Archer | Evaluation of global wind power[END_REF]; [START_REF] Archer | Chap. An Introduction to meteorology for airborne wind energy[END_REF]). De plus, un cerf-volant en vol dynamique aura une vitesse propre, ce qui augmente son vent apparent et donc la force de traction qu'il peut produire.

Pour pouvoir développer de telles ailes, capables de tracter des navires de transport par exemple, des modèles numériques sont nécessaires afin de pouvoir estimer les performances d'un cerf-volant, sa puissance et sa résistance, sans avoir à réaliser et tester différents prototypes, de plus ou moins grande taille. De plus, de nombreux paramètres peuvent influer sur les performances d'un cerf-volant, sa géométrie, les matériaux utilisés pour sa fabrication, le jeu de brides utilisé et leurs points d'accroche sur le cerfvolant, etc, ce qui augmente d'autant le nombre de tests à réaliser. Les cerfs-volants utilisés dans le cadre de ce projet sont des ailes à boudins gonflables (voir Figure 1.2). Ces boudins gonflés permettent de donner une première forme au kite, et de rigidifier légèrement la structure mais le cerf-volant reste une structure souple, dont la déformation dépend du chargement subi. Or, en vol, la vitesse et donc le chargement du kite dépend de sa géométrie. Il est donc nécessaire pour pouvoir estimer correctement les performances de ce genre de cerf-volant de mettre en place des modèles numériques prenant en compte l'interaction fluide-structure.

Aujourd'hui, il est possible de modéliser complètement un cerf-volant en interaction fluide-structure avec un modèle éléments finis pour la structure [START_REF] Leloup | Model approach and numerical tool development for kite performance assesment and mechanical design; application to vessels auxiliary propulsion[END_REF]; [START_REF] Bosch | Dynamic Nonlinear Aeroelastic Model of a Kite for Power Generation[END_REF]) et volumes finis pour la partie fluide [START_REF] Duport | Comparison of 3D non-linear lifting line method calculations with 3D RANSE simulations and application to the prediction of the global load on a cornering kite. 15ème Journées de l'Hydrodynamique[END_REF]; [START_REF] Scupi | Large Eddy Simulation of Flows Around a Kite Used as an Auxiliary Propulsion System[END_REF]), mais ces modèles restent très coûteux en terme de puissance et de temps de calcul, surtout pour mettre en place un couplage en interaction fluide-structure. C'est pourquoi l'objectif de cette thèse est de développer des modèles rapides et fiables, capables d'estimer rapidement les performances d'un cerf-volant afin de pouvoir tester, pour une première phase de conception, différentes géométries d'ailes, dans différentes conditions. Des modèles plus complexes pourront ensuite être utilisé sur quelques cas particuliers, un cas de chargement critique par exemple, afin de finir la conception du cerf-volant.

Le modèle fluide est présenté dans le Chapitre 2. Ce modèle est une extension de la v théorie de la ligne portante de Prandtl, conçu pour pouvoir gérer des ailes arrondies, avec des angles de dièdre, de vrillage ou de flèche variables le long de l'envergure. Le modèle est aussi capable de prendre en compte le taux de rotation du kite, différents angles d'incidence ou de dérapage du kite, ainsi que la non-linéarité du coefficient de portance en se basant sur un algorithme proposé par [START_REF] Anderson | Fundamentals of aerodynamics[END_REF]). La construction des sections de segments tourbillonnaires est basée sur une géométrie proposé par Katz et Plotkin [START_REF] Katz | Low speed aerodynamics[END_REF]). Le tourbillon lié est placé au quart de corde de l'aile et est orthogonal aux deux tourbillons libres adjacents, qui s'étendent parallèle à la corde sur une longueur de corde. Le sillage de l'aile est représenté par deux autres segments tourbillonnaires, parallèles au vent apparent, qui s'étendent sur plusieurs longueurs de corde. Un dernier segment permet de refermer la section (voir Figure 2.5). L'algorithme permet de calculer la circulation de chacun des segments et ainsi la portance, la trainée et le moment, en local le long de l'envergure ou global sur l'aile. Pour la bonne convergence de l'algorithme, un coefficient de relaxation est utilisé. Le programme ne met que quelques secondes pour converger, sans avoir été efficacement optimisé. De nombreux tests numériques permettent de recommander une discrétisation entre 60 et 100 sections sur l'envergure de l'aile (voir Figure 2.11), ainsi qu'une longueur de sillage entre 10 et 20 fois la corde à la racine (voir Figure 2.9), pour une bonne convergence du modèle. Enfin, le modèle coïncide correctement avec les valeurs prévues par la théorie ligne portante sur des géométries d'ailes classiques (aile elliptique avec ou sans vrillage, voir Figures 2. 13 et 2.14).

Ce modèle est vérifié dans le Chapitre 3 grâce à des simulations RANSE 3D, sur différentes géométries d'aile avec des angles d'incidence et de dérapage allant de -4°à 16°e t de 0°à 15°respectivement. Tous les calculs sont faits avec la section NACA2412, différente d'une section aérodynamique de kite à boudins gonflés, mais beaucoup plus facile et moins coûteuse à mailler correctement. Premièrement, des simulations sont faites pour calculer les coefficients aérodynamiques utilisés dans le modèle ligne portante (afin qu'ils correspondent aux simulations RANSE) et pour estimer la précision numérique des résultats produits (voir Figures 3.7,3.8 et 3.6). Les résultats concordent avec moins de 5% d'écart sur l'estimation en portance, pour une aile avec des lois de flèche et de vrillage non nuls, en incidence pure avec des angles compris entre -4°et 16°. L'estimation de la traînée est aussi très correct en gardant à l'esprit la différence de temps de calcul, les résultats diffèrent seulement pour des angles d'incidence supérieur à 10°. Pour l'aile en incidence et en dérapage jusqu'à 15°, le modèle est aussi parfaitement capable de suivre l'évolution des simulations RANSE (voir Figure 3.18). Enfin, ces simulations permettent de vérifier le chargement local calculé par le modèle ligne portante. Le modèle Ligne Portante 3D Non-Linéaire est capable d'estimer correctement la position du minimum et du maximum de chargement, pour des ailes avec des angles d'incidence et de dérapage dans la zone stationnaire (voir Figures 3.19 et 3.20).

Le modèle est ensuite implémenté dans une nouvelle procédure itérative d'équilibrage du kite, qui est utilisée pour estimer le chargement du cerf-volant en vol sur une trajectoire circulaire (voir Figure 4.1). Dans ce cas, l'influence des taux de rotation sur le vi chargement est mise en évidence (voir Figure 4.2). Une étude paramétrique est aussi menée pour évaluer et analyser les variations des coefficients aérodynamiques caractéristiques du vol en fonction de la longueur des lignes et du rayon de la trajectoire (voir Figures 4.3 et 4.4).

Le modèle structure Kite as a Beam est présenté dans le Chapitre 5 ainsi qu'un modèle éléments finis plus complet, le modèle IDK, évolution du modèle proposé par [START_REF] Leloup | Model approach and numerical tool development for kite performance assesment and mechanical design; application to vessels auxiliary propulsion[END_REF]. L'idée centrale de ce modèle est de représenter un cerf-volant à boudins gonflés par un ensemble de poutres (voir Figure 5.5). Le kite est divisé en plusieurs cellules élémentaires, chacune composé d'une partie du bord d'attaque (un boudin gonflé), de deux lattes et de la canopée correspondante. Chacune de ces cellules élémentaires est remplacé par une poutre aux propriétés mécaniques équivalentes. D'un point de vue structurel les boudins gonflés, lattes et bord d'attaque, sont représentés par des éléments poutre dont les propriétés mécaniques sont estimées à partir des propriétés du tissu qui les composent (voir Section 5.1.1). La canopée est elle modélisée par des éléments coque, plus stable numériquement que des éléments membrane (voir Section 5.1.2). Pour le calcul des propriétés mécaniques des poutres équivalentes, la cellule élémentaire est tout d'abord gonflée avec une pression uniforme, ses quatre coins étant encastrés. Différentes perturbations sont ensuite appliquées à la cellule, sous la forme de déplacement élémentaire. Les raideurs équivalentes de la poutre sont calculées à partir des efforts et des moments de la cellule en réaction à ces déplacements élémentaires (voir Section 5.3.2).

Ce modèle structure Kite as a Beam est ensuite comparé à l'autre modèle éléments finis, le modèle IDK. Pour pouvoir comparer correctement ces deux modèles, et comme le modèle Kite as a Beam fonctionne avec une pression uniforme sur l'envergure, le modèle IDK est tout d'abord gonflé avec une pression constante (voir Figure 6.1). Les déformations de la géométrie sont ensuite estimées et appliquées au modèle Kite as a Beam dont les propriétés ont été calculées au préalable avec la même pression. Enfin, les deux modèles sont soumis aux mêmes déplacements élémentaires avec des conditions aux limites supposées équivalentes. Les cas de comparaison se composent de trois déplacements suivant les axes globaux du kite et d'une rotation de l'oreille du cerf-volant, une de ses extrémités (voir Figure 6.3). Les forces et les moments de réaction sont comparés sur deux géométries différentes, l'une étant issu d'un kite conçu par beyond the sea® et l'autre étant beaucoup plus académique, circulaire et à corde constante. Ces deux géométries permettent de tester les différences d'un point de vue structure pour une géométrie de kite classique et aussi d'évaluer l'influence des approximations géométriques faites lors de la construction du modèle Kite as a Beam. Les résultats montrent que le modèle Kite as a Beam surestime grandement la torsion en partie à cause de ses approximations géométriques (voir Figures 6.5 et 6.7). De plus, le modèle Kite as a Beam n'est pas capable de suivre les évolutions en fonction de la pression du modèle IDK pour les déplacements dans le plan orthogonal au plan de symétrie du kite, même si les résultats restent dans le même ordre de grandeur. Enfin, une étude montre l'influence de la géométrie de la cellule élémentaire (le maximum de creux et sa position) sur les propriétés des poutres équivalentes (Figures 6.12 et 6.13). vii Finalement, le couplage des modèles structures et du modèle fluide est présenté dans le Chapitre 7. La modélisation des brides est détaillée dans la Section 7.2.1. L'échange des données (déformations de l'aile et chargement aérodynamique) entre les deux modèles est très simple pour l'interaction fluide-structure avec le modèle Kite as a Beam car les maillages sont équivalents. Pour le modèle IDK, le chargement du modèle structure se fait à partir de l'angle d'incidence effectif calculé le long de l'envergure de l'aile, qui permet de retrouver la distribution de pression équivalente sur le profil, calculée avec XFOIL (voir Section 7.2.2). La comparaison des deux modèles se fait sur un cas relativement simple, où la corde dans le plan de symétrie de l'aile est à incidence fixée. Après convergence de l'itération structure, les déformations de l'aile sont donc calculés en enlevant le mouvement de corps rigide du kite au bout de ses lignes. Le gain en temps de calcul du modèle Kite as Beam est réellement important, le modèle Kite as a Beam ne mettant que quelques minutes à converger alors que le modèle IDK nécessite une ou deux heures. En revanche, la raideur en torsion et la raideur en flexion dans le plan des poutres du modèle Kite as Beam doit être calibré pour pouvoir retrouver les résultats du modèle IDK (voir Figures 7.3 et 7.18). De plus, une étude de convergence en maillage montre que l'interaction fluide-structure avec le modèle Kite as a Beam est sensible à la discrétisation du maillage (voir Figure 7.14). viii 6.5 Gap (in %)

List of Figures

( #» x c , #» y c , #» z c )
Elementary cell coordinate system

(K, #» x b , #» y b , #» z b ) Kite coordinate system (M, #» x s , #» t , #» z s ) Aerodynamic section coordinate system # » F a
Aerodynamic forces on the wing

# » V a
Apparent wind velocity

# » V K Wing velocity #» V e
Effective wind

# » V RW

True wind velocity

A

Leading edge point c

Wing chord length

C D 3D drag coefficient C d 2D drag coefficient C L 3D lift coefficient C l 2D lift coefficient C M 3D moment coefficient C m 2D moment coefficient C p Pressure coefficient c r
Wing chord length at the root Global warming and emissions reduction of greenhouse gas effect are the challenges of our century as greenhouse gas emissions from human activities are the primary cause of global warming [START_REF] Pachauri | Climate Change 2014: Synthesis Report[END_REF]). The conclusion of the special report (IPCC ( 2018)) of the Intergovernmental Panel on Climate Change (IPCC) on the impacts of a global warming of 1.5°C are highly alarming. At the current rate, global warming is likely to reach 1.5°C before mid century and climate-related risk such as drought, floods or extreme storms are projected to be higher for each tenth of degree Celsius. Global warming means long lasting and irreversible changes and even loss of some land ecosystems and thus species extinctions. In the oceans, global warming is projected to induce increases in ocean temperature and acidity, and decreases in ocean oxygen levels which will impact the growth, development and survival of a broad range of species (IPCC ( 2018)). These changes will have an impact on human activities, on economic growth, health, food security or water supply of human communities. A slower and limited global warming enables greater opportunities for adaptation in the human and ecological systems.

To limit global warming to In April 2018, the International Maritime Organization (IMO) has therefore recently adopted a resolution describing an initial greenhouse gas strategy for international shipping. The objectives are the reduction by 40% in carbon intensity by 2030, 70% by 2050 compared to 2008 levels and a full decarbonization by 2100. The measures include the improvement of the existing fleet, the development and the provision of zero-carbon or fossil-free fuels or the reduction of the global ship speed.

Wind as a source of energy

Global warming is forcing us to rethink our consumption of energy, to discover or rediscover energy sources. Wind as a source of energy has been used since at least the 5 th millennium B.C.. Indeed, one of the oldest known representation of a sailing boat is a painted disc found in Kuwait, dating from the 5 th millennium (Figure 1.1) while windmills can be found as early as the first century. Wind is a green energy, readily available. Furthermore, Archer [START_REF] Archer | Evaluation of global wind power[END_REF]) collected data from wind stations around the world and shows the great potential of wind energy if we are able to collect it efficiently. Archer also showed that regions with sufficient wind resources are widespread around the world. In Europe, for example, the regions with the largest potential are the northern coast, near the North Sea or the coast of the United Kingdom. The capacity of offshore wind farm can not be neglected either. Indeed, above the water, the surface roughness is reduced, which increases the wind speed and the wind power production in comparison to onshore plants. Archer points out the fact that offshore stations receive average wind speed which can be 90% greater than land-stations.

The wind power density, the wind power per unit area, depends linearly from the air density and with a cubic dependency on the wind speed [START_REF] Archer | Chap. An Introduction to meteorology for airborne wind energy[END_REF]), thus the importance of the wind speed estimation. Near the Earth's surface (below 500m), to describe analytically the wind speed V Archer [START_REF] Archer | Chap. An Introduction to meteorology for airborne wind energy[END_REF]) used Equation 1.1,

V (z) = V (z ref ) z z ref γ (1.1)
with z the height and γ a friction coefficient varying from 0.1 above calm water, to 0.4 above big cities, z ref a reference height at which a measured velocity V (z ref ) is done. The power law of Equation 1.1 has been proven empirically to be a good approximation for the wind speed below 500m [START_REF] Archer | Chap. An Introduction to meteorology for airborne wind energy[END_REF]).

The wind velocity is therefore faster at high altitude, which means more energy to capture. The wind is also steadier at high altitude and also above water as there is no obstacles. To reduce our fossil fuel consumption, it perfectly makes sense to look again, and with new ideas, at the way of capturing wind energy.

Kites to produce energy?

The first known kites were developed in China around the 5th century B.C. Since then, their use have varied, from military applications to scientific purposes or just as toys. [START_REF] Pocock | A treatise on the aeropleustic art, or Navigation in the air: by means of kites, or buoyant sails: with a description of the charvolant[END_REF]), for example, investigated the potential of kite for towing vehicles as soon as the 19 th century. Hobbs [START_REF] Hobbs | A quantitative study of kite performance in natural wind with application to kite anemometry[END_REF]) summarizes a big part of kite history for scientific or industrial applications between the 18 th century and the end of the 20 th century.

Kites can have many shapes depending on their application (see Figure 1.2). Nowadays, the design of a kite for kite surfing is still mainly done on the basis of trial and error, by looking for a compromise between stability, performance and maneuverability.

As the power potential of the wind increases with the wind speed, it is interesting to look for wind at high altitude, which is steadier and stronger . Kites, contrary to windmills are less limited in altitude because of the tethers length than can be adjusted. For that reason, the idea of capturing wind energy with kites has been investigated since the end of the 20 th century [START_REF] Goela | Performance Characteristics of a Kite-Powered Pump[END_REF]; [START_REF] Loyd | Crosswind kite power (for large-scale wind power production)[END_REF]). Loyd, for example, describes various ways of producing power with kite, a kite in static flight exerting a tensile stress on the tether and unwinding them from a drum, a kite in crosswind motion here again pulling on the tether or by adding an air turbine on the kite in dynamic flight.

The study highlights the power potential of the kite but also points out the engineering challenges for the implementation of such a project. The first idea has been further investigated these last years [START_REF] Fechner | Airborne Wind Energy[END_REF]; [START_REF] Olinger | Performance characteristics of a 1 kW Scale Kite-Powered System[END_REF]; [START_REF] Argatov | Energy Conversion Efficiency of the Pumping Kite Wind Generator[END_REF]), providing lots of developments on kite control or on the modeling of kite dynamics. Another application for towing kites, already mentioned by Pocock in the 19 th century, is the use of kites for towing ships. This idea was popularized during the 1980s by the speed record of the Jacob's Ladder [START_REF] Wellicome | Ship propulsive kites: an initial study[END_REF]), a boat towed by kite; or in the 2000s with the growth of kite surfing and the new speed records.

Overview of the study

The use of kites for towing ships is therefore not a new idea. In the 1980s, Duckworth [START_REF] Duckworth | The application of elevated sail (kites) for fuel saving auxiliary propulsion of commercial vessels[END_REF]) carried out experiments with an 8 ton workboat, using stable parachute sails. Wellicome [START_REF] Wellicome | Ship propulsive kites: an initial study[END_REF]; [START_REF] Wellicome | Some comments on the relative merits of various wind propulsion devices[END_REF]) developed a first model to predict the kite aerodynamic forces and its performance. He concluded that, for a medium-sized commercial vessel (4500 ton ship), the kite needed to ensure a full propulsion of the ship is really large and that the most practical application at the time is for fuel saving wind assistance or emergency use.

Naaijen [START_REF] Naaijen | On the power savings by an auxiliary kite propulsion system[END_REF]) estimated the fuel saving potential of kite towing ships.

As the estimation depends on the wind velocity, he showed the necessity to optimize the shipping route. His estimation of fuel saving can go up to 50% at Beaufort 7 using a kite of 500m 2 attached to a 350m towing line for a 50 000 dwt (dead weight tonnage) tanker. [START_REF] Leloup | A continuous and analytical modeling for kites as auxiliary propulsion devoted to merchant ships, including fuel saving estimation[END_REF]) also implemented a procedure to predict the fuel saving potential by optimizing the elevation of the kite and the trajectory orientation and position (azimuth and elevation). The prediction are higher than Naaijen's because Leloup introduced additionally the modeling of static flight in the optimization process, when Naaijen only took into account the dynamic flight. The kite can also switch between horizontal and vertical flight paths, which is useful for upwind conditions. Leloup predicted a fuel saving of about 10% for a 50 000 dwt tanker using a kite of 320m 2 with a wind velocity of 10m.s -1 .

The German company SkySails was one of the first to equip a cargo ship with a kite as an auxiliary propulsion device. In 2008, the company set up a prototype of kite propulsion system of 320m 2 on a 132m ship [START_REF] Erhard | Control of Towing Kite for Seagoing Vessels[END_REF]; [START_REF] Fritz | Airborne wind energy[END_REF]).

In 2017, the company equipped the 35m catamaran Race For Water with a 40m 2 ram air kite (see Figure 1.3). An other vessel towed by kite was launched in 2017, Energy Observer (see Figure 1.3), but this time equipped with a Leading Edge Inflatable (LEI) kite from the beyond the sea® project.

Presentation of the beyond the sea® project

The beyond the sea® project started in 2007 with the aim to develop kite as an auxiliary propulsion for any kind of vessels, to reduce the environmental impact of shipping. The advantages of the kite for the use of wind energy are numerous. The kite and its fixation take less room on the deck than a rig, the heeling moment of a kite is much lower than a sail. Kite can also fly higher, where the winds are stronger and steadier, thanks to adjustable tether lengths. Furthermore, kites have their own velocity in dynamic flight, which increases the apparent wind, and then the towing forces. The inflatable leading edge is an other advantage, it shapes the kite and makes the launch easier.

The project got the support of the French Environment and Energy Management Agency (ADEME) in 2014 and formed a consortium made of the following companies: beyond the sea® , Porcher Industries, Cousin Trestec, BOPP, CMA-CGM and ENSTA Bretagne. The beyond the sea® company manages the various aspects of kite design and automatic control, as well as the launching and recovery procedures. This problematic led to the thesis of Du Pontavice (Du Pontavice (2016)) for example. The technical fabrics needed for the kite conception are made by Porcher Industries. Indeed, fabrics as light as possible are needed, so that the kite can fly with a minimum wind, but it must also resist to the large aerodynamic loadings, as well as to maritime environment exposure. Cousin Trestec manufactures the tethers and bridles used in the beyond the sea® project. Here again, it needs to be as light as possible but also able to resist to the maritime environment. BOPP provides the deck equipment needed to control the kite and to fix it to the deck. CMA-CGM, the third biggest shipowner in the world, can provide appropriate ships for full scale tests. Eventually, ENSTA Bretagne is a research partner of beyond the sea® since early 2010 and the PhD program of Richard [START_REF] Leloup | Model approach and numerical tool development for kite performance assesment and mechanical design; application to vessels auxiliary propulsion[END_REF]).

This first thesis [START_REF] Leloup | Model approach and numerical tool development for kite performance assesment and mechanical design; application to vessels auxiliary propulsion[END_REF]) has investigated various research axes, such as the estimation of the fuel saving potential or a first fluid-structure model of the kite. It was followed by two other PhD programs, expanding further some problematics. Nedeleg Bigi [START_REF] Bigi | Investigation of the dynamic motions and operability of a ship towed by kite[END_REF]) focused on the interactions between ship and kite. Indeed, the kite can impact the maneuverability of the ship and the waves can impact the stability of the kite. For this reason, a dynamic model of the system formed by the kite and the ship is necessary. Bigi showed the importance of taking into account the motion of the ship and its influence on the kite flight.

The second PhD program was pursued by Morgann Behrel [START_REF] Behrel | Investigation of kites for auxiliary ship propulsion: experiment setup, trials, data analysis and kite specs novel identification approach[END_REF]). The aim of his work was to produce experimental data to check the numerical models. A test bed has been developed which allows the control of the kite and the measures of various parameters, such as the tension in the tethers. Several measurement campaigns have been done, including one on land with the test bed fixed on the ground and one at sea, on a boat specifically build for this purpose.

Introduction to the thesis

The beyond the sea project® thus aims to equip ships with kites of various shapes and sizes. Kites far larger than the ones used for common kite surfing need to be designed. For such large kites, a design method based of trial and error with real testing is not possible anymore. Therefore a need arose for numerical models able to determine the kite performances and the loads in the structure. As a kite is a flexible structure, these models have to take into account the fluid-structure interactions (FSI). Indeed, the aerodynamic performances of a kite depend on its geometry, which depends on its loading. Coupling the fluid and the structure problems is a necessity to get the flying shape of the kite.

Even if the structural problem can be modeled by a complete Finite Element (FE) model [START_REF] Leloup | Model approach and numerical tool development for kite performance assesment and mechanical design; application to vessels auxiliary propulsion[END_REF]) and the fluid problem by a Computational Fluid Dynamics (CFD) method, a fully coupled FSI simulation using FE and CFD methods is still highly computationally demanding. It can then be very useful in a design stage to use a fast, reliable and simple model to estimate rapidly the kite aerodynamic performances. Moreover fast simulations allow the studies of different wind conditions and the determination of the critical loading case, where more complex and more time-consuming models can't be used so easily.

A fast FSI model should be able to compute the geometrical deformations of the wing, the tensions in the tethers, the aerodynamic loading distribution and the aerodynamic coefficients of the wing by taking into account the geometry of the wing, the characteristics of the wing fabrics; the mechanical properties of the tethers or various sets of bridles. As the kite can follow a complex trajectory, the method has also to take into account the wind direction, the wing angle of incidence or its sideslip angle, or the turn rate of the kite. A fast model can be used in a first design stage to test rapidly different kite geometries, shapes and sizes, or different sets of bridles as well as different wind conditions. In this thesis, it was chosen to build the fluid model and the structure model separately. In that way, each model can be used or improved independently.

The present document is divided into three main parts. The first part (Chapters 2, 3 and 4) focuses on the fluid model. Chapter 2 presents the fluid model algorithm and some numerical tests, such as an estimation of the numerical accuracy. Chapter 3 is devoted to the comparison of the fast fluid model with RANSE simulations, in order to validate the obtained results. Chapter 4 presents an use of the fluid model to simulate the kite in flight as a rigid body.

The second part (Chapters 5 and 6) deals with the structural problem. Two structure models are presented in Chapter 5, the fast one developed over this thesis and an evolution of a complete Finite Element model developed in the beyond the sea® project. Chapter 6 describes some elementary comparison cases carried out to check the differences between the two models.

Finally, the last part and last Chapter of this document is the implementation of the fluid-structure interactions. Chapter 7 presents the coupling between both models and few computation cases.

2 Fluid model presentation

Introduction to the Lifting Line problem

The 3D Non-Linear Lifting Line model is based on the classical Prandtl's lifting line theory, which is summarized in this Section. For more details, refer to Katz and Plotkin [START_REF] Katz | Low speed aerodynamics[END_REF]), [START_REF] Anderson | Fundamentals of aerodynamics[END_REF] or [START_REF] Bertin | Aerodynamic for engineers[END_REF]).

The Prandtl's lifting line theory models a wing by a number of horseshoe vortices. Each horseshoe vortex is composed of one bound vortex placed at the quarter chord of the wing section along the span. Since a vortex filament cannot end in a fluid, according to the Helmholtz's vortex theorem, two free vortices aligned with the face stream finish the horseshoe vortex (see Figure 2.2). The strength of the vortex filament is defined by its circulation Γ which is constant all along, in accordance with the second part of the Helmholtz's vortex theorem.

The velocity induced by the vortex filament is given by the Biot-Savart law (see [START_REF] Karamcheti | Principles of Ideal-Fluid Aerodynamics[END_REF]). For example, Equation 2.1 gives the velocity induced at point M by the vortex filament (L) of Figure 2.1 with P a point on (L), # » dL P = dL P #» t P and #» t P the unit tangent vector at point P .

# » V ind = Γ 4π (L) #» t P × # » P M # » P M 3 (2.1)
For a planar wing with a planar wake (case presented in Figure 2.2), the induced velocity at point M , # »

V ind (M ) is given by computing Equation 2.1 for each vortex segment.

# » V ind (M ) = - 1 4π b/2 -b/2 Γ (y) y 0 -y dy #» z (2.2)
The Kutta-Joukowski theorem states that lift per unit span on a two-dimensional wing is directly proportional to the circulation around the wing:

L = ρ # » V ∞ × Γ # » dy = ρV ∞ Γdy (2.3)
Therefore, in Prandtl's lifting line theory, the determination of the circulation provides the aerodynamic specs of the wing, the lift and induced drag, and also the induced wind on the wing, which gives its effective angle of incidence.

2 -b 2 ∞ ∞ ∞ ∞ ∞ Γ(y + dy) Γ(y) #» x #» y #» z M (0, y 0 , 0) P (0, y, 0) Γ(y + dy) -Γ(y) = Γ (y)dy
Lifting Line model Graf [START_REF] Graf | Comparison of full 3D-RANS simulations with 2D-RANS/lifting line method calculations for the flow analysis of rigid wings for high performance multihulls[END_REF]) uses a non-linear iterative lifting line method similar to the method described by [START_REF] Anderson | Fundamentals of aerodynamics[END_REF]) to predict the lift and drag of a twoelement straight wing for AC72 catamaran. This method is not based on the linear relationship between the section lift coefficient and the section angle of attack and for that reason it could be used to take into account the effects of stall. The comparison with RANSE simulations shows a good agreement for attached flow regime and in the first part of the nonlinear range before the maximum lift. Furthermore, the method is computationally efficient enough for some optimization. However, in case of larger angles of incidence, the lifting line method shows its limits as it is sometimes unable to converge. This may be a problem for the modeling of the kite in extreme flight conditions.

Phillips [START_REF] Phillips | Modern adaptation of Prandtl's classic lifting-line theory[END_REF]) also developed a non-linear lifting line method able to deal with wings with both sweep and dihedral angles. The method is validated for classic elliptic wings and then compared with results obtained from a numerical panel method, from an inviscid CFD solution and from experimental data. The results shows a good agreement between the CFD and the lifting line method, which shows the possibilities of the lifting line method for wings with dihedral and sweep angles.

Gaunaa [START_REF] Gaunaa | A computationally efficient method for determining the aerodynamic perfomace of kites for wind energy applications[END_REF]) developed a computationally efficient method to determine the aerodynamic performances of kites. The approach iteratively couples a Vortex Lattic Method to 2D airfoil data to consider the effects of airfoil thickness and of viscosity. The results are compared with RANSE simulations and show a good agreement for cases only in incidence and without stall. Gaunaa also presents results for a kite in translation with a sideslip angle, which are consistent, at least for low angles of sideslip.

Leloup [START_REF] Leloup | Estimation of the Lift-to-Drag ratio using the lifting line method: application to a leading edge inflatable kite[END_REF]; [START_REF] Leloup | Model approach and numerical tool development for kite performance assesment and mechanical design; application to vessels auxiliary propulsion[END_REF]; [START_REF] Leloup | Estimation par la méthode de ligne portante de l'effet d'un virage sur la finesse: application aux cerfs-volants de traction pour la propulsion auxiliaire des navires[END_REF]) adapted the lifting line method for 3D kite wings with variable dihedral and sweep angles along the span. However this strategy is a direct solving based on a collocation method, which prohibits to take into account the non-linearity of the lift coefficient. This limits the use of this lifting line method to small angles of incidence.

This Chapter describes the 3D non-linear lifting line model developed to enable the estimation of the aerodynamic performances of a wing with variable dihedral and sweep angles. Section 2.2 explains the construction of the lifting line, the algorithm used for the circulation computation and the post-processing, the determination of the aerodynamic coefficients. Section 2.3 presents the influence of the initial parameters, the discretization or the wake length. An estimation of the numerical accuracy of the method is presented, as well as a first verification of the results thanks to the analytical solution of the elliptic wing.

3D Non Linear Lifting Line model

This part is devoted to the presentation of the 3D non-linear lifting line algorithm. Section 2.2.1 is committed to the presentation of the lifting line geometry while Section 2.2.2 presents the way chosen to solve the circulation distribution. Section 2.2.3 details the computation of the aerodynamic coefficients of the wing.

Parametric definition of the lifting line

The kite geometry is defined in the Cartesian reference frame

(K, #» x b , #» y b , #» z b ) (see Figure 2.3). A generatrix line is given as a parametric planar curve M(s) in the plane (K, #» y b , #» z b ),
M denoting the current point and s the curvilinear abscissa. This curve is elliptical and symmetrical with respect to the axis (K, #» z b ). At the current point M(s), the section of the kite is defined in the normal plane with respect to the tangent vector to the generatrix line

#» t(s) = # » dM /ds. This plane is denoted by (M(s), #» x b , #» n(s)), where #» n(s) = #» x b × #» t(s)
is the normal vector to the generatrix line in the (K, #» y b , #» z b ) plane. The local chord direction of the section #» x s (s) is obtained from the rotation of the vector #» x b of an angle α v (s) around #» t(s), where α v (s) is a given twist law. Since M(s) is supposed to be the quarter chord point of the current section, the leading edge A(s) and the trailing edge F(s) are then located along the chord axis (M(s), #» x s (s)) according to a given chord law c(s), following s). This finally leads to the local chord reference frame of the kite section (A(s), #»

# » AM(s) = 0.25c(s) #» x s (s) and # » M F(s) = 0.75c(s) #» x s (
x s (s), #» z s (s)) with #» z s (s) = #» x s (s) × #» t(s),
where the points of the extrados and intrados can be placed; for example from a given non dimensional section definition, scaled with the chord law c(s). In the case of a kite wing with sweep angle, the local section points are finally translated by the vector

f(s) #»
x b , with f(s) a given sweep law. Appendix 7.4.3 details the different twist, sweep and chord laws implemented in the house made LP3DNL program as well as the various geometrical parameters used for the construction of the geometry.

#» n #» zs αv #» y b #» t M #» x b # » xs αv K #» z b #» x b M c # » xs f #» x b αv A F M Figure 2.

-Definition of the lifting line

Lastly each section must be coupled to its aerodynamic properties. For the model, these aerodynamic specs consist of three 2D polar curves, lift, additional drag and moment, with respect to the angle of incidence. These 2D polar curves are described by a combination of polynomials of various degrees, which have been chosen in order to make easier the copy of polar curves obtained via experimental measures or CFD simulations.

-1 0 1 2 3 (α 1 , C l1 ) (α 2 , C l2 ) (α 3 , C l3 ) (α 4 , C l4 ) I II III IV V Angle of incidence (°) Lift Coefficient

The lifting line method

The 3D non-linear lifting line is based on an extension of the Prandtl's lifting line theory. This extension is intended to address cases of wings with variable dihedral and sweep angles. Leloup introduced this technique in [START_REF] Leloup | Estimation of the Lift-to-Drag ratio using the lifting line method: application to a leading edge inflatable kite[END_REF] with a linear implementation while the present method is taking into account the non-linearity of the lift coefficient, generalizing the work of [START_REF] Anderson | Fundamentals of aerodynamics[END_REF]) and Graf [START_REF] Graf | Comparison of full 3D-RANS simulations with 2D-RANS/lifting line method calculations for the flow analysis of rigid wings for high performance multihulls[END_REF]) for straight wings. Assuming a rigid kite, the wing is supposed to fly in a given wind # »

V RW , with a given velocity # »

V K at the quarter chord point K in its symmetry plane, and with a given turn rate #» Ω. In this general case, the kite global apparent wind is

# » V a = # » V RW - # »
V K and its direction is denoted by # » x a . The finite wing and its wake are represented by a set of horseshoe vortices of different strengths Γ = (Γ i ) i=1...2n . The aim of the algorithm presented thereafter is to calculate the circulation Γ i of each horseshoe vortex. Once these strengths are obtained, the local effective flow for each wing section allows local aerodynamic forces and torques calculation along the span of the wing. The numerical iterative solution is taken from Anderson [START_REF] Anderson | Fundamentals of aerodynamics[END_REF]), but the calculation of the local effective angles of incidence is adapted to the cases of wings which are non-straight and non-planar. The horseshoe vortices used for discretization, and the calculation of their influences, are for their part derived from [START_REF] Katz | Low speed aerodynamics[END_REF]).

The wing is divided in a finite number 2n of plane sections normal to #» t i , as defined in section 2.2.1, each one is represented by a horseshoe vortex. The horseshoe vortex number i consists of six vortex segments. The bound vortex Γ i #» t i l i is located at the local quarter chord length, perpendicular to the plane of the considered section. Each of the two trailing vortices are separated into two parts: the first which extends parallel to the chord over one chord length, ±Γ i #»

x s i c i , and the second which extends parallel to the global relative free stream direction # » x a and over several chords length leading to L wake length, ±Γ i # »

#» V RW #» V K #» Ω K #» x b #» z b #» y b
x a L wake . The starting vortex -Γ i #» t i l i is finally used to close the horseshoe vortex. An example of the discretized model is presented in Figure 2.5 and the geometry of the section is described in Figure 2.6. It can be noticed that even with a swept wing, the bound vortex along the lifting line is orthogonal to the two adjacent trailing vortices. This leads to a piecewise constant discretization of the lifting line as it can be seen in Figure 2.5, but a correct match is needed between the local lift calculated from the Kutta formula and from the polar of the section, as theoretically required by such models.

M A F P 1 P 2 P 3 P 4 P 5 P 6 l i #» t #» x s c #» x a L wake Figure 2.6 -Section geometry
For the section i of the kite, the local apparent wind is given by Equation 2.4 (see Figure 2.5 for the notations).

# » V a i = # » V RW -( # » V K + #» Ω × # » KM i ) (2.4)
The numerical iterative solution for the circulation Γ i of each horseshoe vortex is described hereafter with p the counter of the recursive steps. The local circulation values are first initialized by an elliptical or other distribution along the wing span (see Section 2.3.1), denoted Γ i (0) for p = 0. At step p, starting from the

Γ i (p) values, the Γ i (p+1)
values are updated as follows.

For each point M i of the lifting line, the induced velocities by each vortex segment of the whole wing model are calculated from Γ i (p) values with the Biot-Savart law and then summed, leading to the induced velocity

# » V ind i (p) . # » V ind i (p) = 2n j=1 Γ j (p) 1 4π #» v ij 53 + #» v ij 31 + #» v ij 12 + #» v ij 24 + #» v ij 46 + #» v ij 65 (2.5) with #» v ij kl = # » P j k P j l • # » P j k M i # » P j k M i - # » P j k P j l • # » P j l M i # » P j l M i # » P j k M i × # » P j l M i # » P j k M i × # » P j l M i 2
In Equation 2.5, it can be noticed that only the circulation depends of the step p. The value of the second part of the equation is thus computed only once, at the beginning of the algorithm, and then stored in order to save some computation time.

The local effective wind is obtained by combining the induced wind and the local free stream velocity:

#» V e i (p) = # » V a i + # » V ind i (p) (2.6)
This velocity is projected in the plane of the section, leading to

# » V prj i (p)
and to the section angle of effective incidence

α i prj (p) (see Figure 2.7). # » V prj i (p) = #» V e i (p) -( #» V e i (p) • #» t i ) #» t i (2.7)
Using the 2D lift polar curve of the section, one obtains the local lift per unit length dL i (p) , in the section plane, orthogonal to

# » V prj i (p) . dL i (p) = 0.5ρ(V i prj ) 2 (p) c i C i l (α i prj (p) )
(2.8)

By construction this local lift is orthogonal to the (

#» V e i (p) , #» t i
) plane, as it should be according to the Kutta formula

# » dL i (p) = ρ #» V e i (p) × Γ i (p) #» t i = ρ # » V prj i (p) × Γ i (p) #» t i applied to the bound vortex.
The value of the bound vortex strength Γ i (p+1) is finally settled so that the Kutta formula also leads to the previous modulus dL i (p) . In other words, Γ i (p+1) is calculated by solving Equation 2.9 for the unknown

Γ i (p+1) . ρ # » V prj i (p) × Γ i (p+1) #» t i = dL i (p)
(2.9) Equation 2.9 naturally leads to:

Γ i (p+1) = 0.5V i prj (p) c i C i l (α i prj (p) )
(2.10)

Lifting line M # » V e ( # » V e • #» t ) #» t #» V prj #» t # » dL #» z s #» x prj α prj α v #» x b #» x s Γ Figure 2.7 -Local effective wind projected in the section plane
The use of a damping factor is necessary for the convergence of the algorithm. Anderson [START_REF] Anderson | Fundamentals of aerodynamics[END_REF]) recommends for example a damping factor in the order of magnitude of 0.05. In the present lifting line algorithm, this factor is self adaptive in order to reduce the computational time. If the difference between Γ i (p+1) and Γ i (p) decreases, the damping factor is increased and if this difference increases, the damping factor is decreased following Equation 2.11. with w the damping factor, w min and w max the interval of the damping factor and e w the coefficient used to vary the damping factor. For more details on the adjustment of the user defined values w min , w max and e w , refer to Section 2.3.1.

if 2n i=1 (Γ i (p+1) -Γ i (p) ) 2 2n < 2n i=1 (Γ i (p) -Γ i (p-1) ) 2
The circulation value is ultimately updated by weighting between current Γ (p+1) and previous Γ (p) values using the damping factor:

Γ (p+1) ← Γ (p) + ω Γ (p+1) -Γ (p)
(2.12)

This whole process is repeated and, as expected, the circulation distribution is found to converge according to a stopping criteria of the form:

2n i=1 (Γ i (p+1) -Γ i (p) ) 2 2n <   10 * 2n i=1 (Γ i (p+1) ) 2 2n + 1   * 10 -(ncs+1) (2.13)
where n cs is approximately the number of converged digits required by the user.

Post-processing

After the convergence of the algorithm, the circulation is known on each section of the wing. The next step is to calculate the local loads which are needed for the fluidstructure interaction. For the computation of the local lift effects # » dR l i , two methods were implemented, one using the Kutta formula (Equation 2.14) and one using the lift polar of the section (Equation 2.15).

The Kutta formula gives directly the load in the kite coordinate system with #» V e i the effective wind on the section, Γ i its converged circulation and l i the width of the section.

# » dR l i = ρΓ i #» V e i × l i #» t i (2.14)
For the calculation of the local lift with the lift polar curve, the effective wind is projected in the section plane in order to get the meaningful effective wind V i prj , and the effective angle of incidence of the section α i prj . By definition, # » dR l i is perpendicular to the projected effective wind.

dR i l = 1 2 ρ V i prj 2 c i C i l (α i prj )l i (2.15)
Similarly, the local additional drag effects are computed with the drag polar curve, Equation 2.16 and they are collinear to the projected effective wind.

dR i d = 1 2 ρ V i prj 2 c i C i d (α i prj )l i (2.16)
The local moment at the quarter chord is orthogonal to the plane of the section and its magnitude is given by Equation 2.17:

dM i = 1 2 ρ V i prj 2 c i 2 C i m (α i prj )l i
(2.17)

The global aerodynamic loads are then calculated by summing the local loads :

# » F a = 2n i # » dR l i + # » dR d i (2.18)
The global drag

#» D is collinear to the apparent wind # » V a while the lift #» L is orthogonal to the apparent wind (Equation 2.19). #» D = # » F a • # » x a # »
x a and

#» L = # » F a - #» D (2.19)
The global moment of the wing expressed at point K, point at the quarter chord in the middle of the wing is given by Equation 2.20:

# » M K = 2n i # » dM i + # » KM i × # » dR l i + # » dR d i (2.20)
The aerodynamic coefficients of the wing can now be computed with S the surface of the wing and c r the chord length at the root of the wing, in the symmetry plane:

C L = || #» L|| 1 2 ρSV 2 a (2.21) C D = || #» D|| 1 2 ρSV 2 a (2.22) C M = || # » M K || 1 2 ρc r SV 2 a (2.23)
The lift-to-drag ratio is conventionally defined by Equation 2.24.

ε = L D = C L C D (2.24)

Numerical tests

This part consists of a series of numerical tests and checks of the 3D non-linear lifting line program. Section 2.3.1 describes the influence of the initial parameters and of some computation parameters such as: initial circulation, damping factor, wake length and lift per section post processing. Section 2.3.2 gives an estimation of the numerical accuracy of the method and Section 2.3.3 compares the results with the classical and analytical Prandtl's lifting line theory on an elliptic wing.

Initial conditions and computation parameters Initial circulation distribution

The algorithm needs to be initialized with a circulation distribution along the wing span. Three different initial circulation distributions have been tested, see Equations 2.25, 2.26, 2.27, Figure 2.8 and Table 2.1. Equation 2.25 is the circulation distribution derived from the Kutta-Joukowski theorem applied to each section and neglecting the 3D effects. Equation 2.26 is the circulation distribution derived from the circulation at the root, weighted by an elliptic distribution along the span. Equation 2.27 is a circulation distribution constant along the span, computed as the average of the local circulation values from Equation 2.25.

the average of the local circulations, constant along the span.

Γ i (0) = 1 2 V i prj c i C i l (α i prj )
(2.25)

Γ i (0) = 1 2 V n prj c r C n l (α n prj )g coef f 1 - (s i -s max ) 2 s max 2
(2.26)

Γ i (0) = 1 2n 2n k=1 1 2 V k prj c k C k l (α k prj k ) (2.27) With V i prj the apparent wind # » V a
i projected in the section plan, c i the section chord, s i its curvilinear abscissa, α i prj the effectif angle of incidence of the section (which depends of the twist and of the apparent wind of the section), n the number of sections on the half-wing and g coef f a user defined coefficient (often g coef f = 0.5). It is found that all these initial circulation distributions have no influence on the converged circulation distribution, the differences between each case bring of the same order of magnitude than the numerical error, see Figure 2.8 where no notable differences can be observed between the converged curves. However the initial circulation distribution has an influence on the number of iterations needed to obtain the converged solution, as shown in Table 2.1. The fastest distribution is found to be the closest of the converged solution. Therefore the initial elliptic distribution of circulation (Equation 2.26) with a coefficient g coef f = 0.5 is recommended .

Kutta-Joukowski theorem (1) Weighted elliptic distribution (2) Mean value (3)
Iterations number 147 135 137

Adaptive damping factor

As explained in Section 2.2.2, an adaptive damping coefficient is used to achieve the convergence of the iterative calculation of the circulation distribution. It is found that the initial damping coefficient has no influence on the convergence outside of the number of iterations if chosen very small (< 0.1). The limits w min and w max of the coefficient and the parameter e w are by contrast really significant. w max can be close to 0.9 but w min has to be lower or equal to 10 -2 (ideally 10 -3 ) for the correct convergence of the algorithm. Similarly, a parameter e w of 0.1 has proved to be very robust and effective during the different uses of the algorithm. Whatever the values chosen for the relaxation bounds and factor w min , w max and e w , they have no influence on final results if the algorithm is correctly converged.

Lift per section post processing

As stated in Section 2.2.2 the geometry of the lifting line model is made in order to have a correct match between the local lift calculated from the Kutta formula and from the polar of the section. This point was checked, and for any sideslip or incidence angle, the difference between the two calculation methods is in the order of magnitude of the numerical error (10 -10 ).

Wake length

To evaluate the influence of the wake length, the aerodynamic coefficients of a Leading Edge Inflatable wing have been computed with a wake length ranging from c r to 1000c r .

The results are shown in Figure 2.9.

10 0 10 1 10 2 10 3 10 -8 10 -6 10 -4 10 -2 N = Wake length/Chord Relative error (%) Lift Coefficient Drag Coefficient Moment Coefficient Figure 2.9 -|C X (N )-C X (10 3 )| |C X (10 3) |
Relative error on the aerodynamic coefficients with respect to the wake length nondimensionalized by the root chord length.

For L wake ∈ [10c; 20c], the three aerodynamic coefficients have at least 3 continuous shapes which is satisfactory taking into account the intrinsic accuracy of the lifting line model. It can be noticed that the wake length has almost no influence on the number of iterations.

Numerical accuracy estimation

Eventually, the influence of the mesh has to be evaluated. First, it is necessary to check the influence of the sections number per half-wing on the computation time. The lifting line program has been run for a wing with between 10 and 500 sections, 50 times each. Figure 2.10 (a) shows the averaged computation time for 50 runs on a classic computer. The purpose here is not to give a precise estimation of the computation time as the program is not properly optimized but more to give an idea of the evolution of the computation time as a function of the number of sections on the wing. The number of iterations can be seen in Figure 2.10 (b).

As can be seen in Figure 2.10, the computation time and the number of iterations start to increase significantly for more than 100 sections per half-wing. Under 100 sections, the influence of the iterations number is harder to evaluate as the other parts of the program, especially files writing, take comparatively more time. Furthermore, as shown in Figure 2.11, for 30 sections per half-wing the estimation of the aerodynamic coefficients has three significant figure which is considered enough, taking into account the intrinsic accuracy of the lifting line model. Indeed, for a large number of sections, the numerical accuracy is certainly superior to the accuracy of the lifting line model itself. That is why, for the use of the lifting line program, a discretization between 30 and 50 sections per half-wing is recommended.

10 1 10 2 10 -5 10 -4 10 -3
Sections number

Relative variation

Lift Coefficient Drag Coefficient Moment Coefficient

Figure 2.11 -Evolution of the relative variation of the aerodynamic coefficients with respect to the number of sections per half-wing. The relative variation is estimated as

|C X (N )-C X (N -1)| 10|C X (N )|+1
, using the results of two consecutive cases, N -1 and N .

For a number of sections per half-wing in the range of [30; 100], the standard deviation relative to the mean value is estimated, which gives the 95% confidence interval. This interval is evaluated for the linear range (angle of incidence typically lower than 10°) at 2°of incidence and for the non-linear range (angle of incidence greater than 10°) at 12°o f incidence. It is assumed that results from these two angles are sufficient to estimate the numerical precision of the lifting line model. In the same way, results obtained with this geometry and this aerodynamic section are resumed in

Verification again an analytical case

The first verification case is the elliptic wing of Prandtl (see Figure 2.12). This wing has an elliptic circulation distribution of the form:

Γ(y) = Γ max 1 - y d 2 2
(2.28)

#» z b #» y b #» x b d 2

Circulation

Leading edge

Trailing edge

Horsheshoe vortices

Figure 2.12 -Wing with an elliptic plan form and an elliptic circulation distribution

For more details on the calculation of the analytical solution, refers to Katz and Plotkin [START_REF] Katz | Low speed aerodynamics[END_REF]) or Anderson [START_REF] Anderson | Fundamentals of aerodynamics[END_REF]). The lift of the wing is given by Equation 2.29 [START_REF] Anderson | Fundamentals of aerodynamics[END_REF]).

L = ρV a d 2 -d 2 Γ(y)dy = πρd 2 V a Γ max 2 (2.29)
Which gives the lift coefficient:

C L = L 1 2 ρSV 2 a = πd 2 Γ max SV a
(2.30)

Katz and Plotkin [START_REF] Katz | Low speed aerodynamics[END_REF]) give the expression for Γ max :

Γ max = 4d 2 V a (α i -α L0 ) 1 + 8d 2 m 0 cr (2.31)
with α i the angle of incidence of the wing, α L0 the zero-lift angle and m 0 the lift coefficient slope (α i belongs to the linear part of the section lift coefficient).

- For the first validation test, the wing has a root chord c r = 1m and a span of 8m, d 2 = 4m, the section is assumed constant along the span. The angle of incidence varies between -5°and 10°, which corresponds to the linear part of the lift coefficient. Figure 2.13 (a) shows the circulation distribution along the span of the wing and Figure 2.13 (b) presents the wing lift coefficient as a function of the angle of incidence for the wing with 30 sections per half-wing (n = 30). The results match admirably, for 10°o f incidence, the gap between the analytical solution and the 3D non-linear lifting line method is of 0.5%, which is within the margin of error of the method.

The second validation case is still an elliptical wing but with a linear twist [START_REF] Katz | Low speed aerodynamics[END_REF], [START_REF] Filotas | Solution of the lifting line equation for twisted elliptic wings[END_REF]). For this case, the circulation is described by the following Fourier expansion (Equation 2.32) with y = d 2 cos(θ), θ ∈ [0; π], which gives the linear twist

α v = α 0 |y/d 2 | = α 0 | cos(θ)| and the elliptic chord c(θ) = c r sin(θ). Γ(θ) = 4d 2 V a ∞ k=1 A k sin(kθ) (2.32)
Using this circulation expression, the equation for the coefficients A k using lifting line theory is:

∞ k=1 A k sin(kθ) 8d 2 m 0 c r sin(θ) + k sin(θ) = α i -α L0 + α 0 | cos(θ)| (2.33) ∞ k=1 A k sin(kθ) 8d 2 m 0 c r + k = (α i -α L0 + α 0 | cos(θ)|) sin(θ) (2.34)
With the Fourier series principle, one can express the coefficients A k from Equation 2.34:

A k = 2 π m 0 c r 8d 2 + km 0 c r π 0 (α i -α L0 + α 0 | cos(θ)|) sin(θ) sin(kθ)dθ (2.35)
Katz and Plotkin [START_REF] Katz | Low speed aerodynamics[END_REF]) demonstrate that the wing lift coefficient depends only of the first term of the Fourier expansion:

C L = 4πd 2 2 A 1 S = 16d 2 m 0 π(8d 2 + m 0 c r ) π 0 (α i -α L0 + α 0 | cos(θ)|) sin 2 (θ)dθ
(2.36)

C L = 16d 2 m 0 π(8d 2 + m 0 c r ) π 0 (α i -α L0 ) sin 2 (θ)dθ + 2 π 2 0 α 0 cos(θ) sin 2 (θ)dθ C L = 16d 2 m 0 π(8d 2 + m 0 c r ) (α i -α L0 ) π 2 + 2α 0 3
(2.37)

Figure 2.14 shows the spanwise circulation (a) and the lift coefficient (b) of an elliptic wing with a maximal twist of 5°in red and -5°in blue. The angle of incidence varies between [-5°; 5°] and [0°; 10°] in order to stay in the linear range of the section lift coefficient. This is enough to confirm the good agreement between the 3D non-linear lifting line method and the lifting line theory for these classic wing geometries.

Fluid model comparison against RANSE simulations

In order to validate the 3D non-linear lifting line method, results were compared to those obtained from 3D RANSE simulations. All the necessary simulations were performed with the CFD code STAR-CCM+®. As a first step, 2D simulations of a NACA2412 wing section were used to adjust the settings of the simulations and for convergence studies, but also to obtain the 2D polar curves of the section, which are required to feed the 3D non-linear lifting line model. In a second step, 3D simulations of various kite geometries were used to calculate local and global aerodynamic forces on a kite in translation, including variations of the angle of incidence and the sideslip angle.

Numerical settings

Within this study, all RANSE simulations were conducted in an incompressible, steady and fully turbulent framework. The retained turbulence model was the two-equation krealizable model with a two-layer formulation for the wall treatment. The segregated flow solver was based on the SIMPLE algorithm, and a second-order discretization scheme. It was chosen to work with a kite section or a kite wing fixed with respect to the computational domain, and to enable the direction of the inlet velocity to vary in order to model the angle of incidence and the sideslip angle. The inlet boundary conditions were such that the chord based Reynolds number was 3.1 * 10 6 , for the root section of the kite. The turbulence intensity was set to 0.5% and the turbulent viscosity ratio was set to 1.

The meshed computational domain is presented in Figures 3.1 and 3.2 for the 2D simulations and in Figure 3.3 for the 3D simulations. The domain was meshed using the trimmed cell mesher, which led to predominantly hexahedral mesh (Figure 3.1). It was controlled by a cell base size (0.1c r ), and targeted cell sizes at some boundaries: inlet and outlet (0.5c r ), wing extrados and intrados (0.025c r ), wing leading and trailing edges (0.00625c r ). Cells size growth rate was very slow, which means at least 8 cell layers of same size per transition (Figure 3.2 (a)). Around the wing a prism layer mesh was used in order to get orthogonal cells next to the wall (Figure 3.2 (b) and (c)). It was controlled by its thickness (0.0125c r ), a number of layers (10), and a growth rate between adjacent cells in the wall normal direction (1.2). For the 2D simulations two For 3D simulations, wake refinements were prescribed on the cylinder which is supported by the trailing edge and whose axis is parallel to the mean free stream, over the same extent as in the 2D case. In order to partially resolve the tip vortex, a supplementary refinement was also prescribed in a conical region, which includes the wing tip and whose axis is parallel to the mean free stream. These refinements are illustrated The tolerance was set to 10 -6 for both coefficients, and it was found that it corresponded to the drop of the non-dimensional residuals over at least 4 or 5 decades.

These simulations were also necessary to check the local results of the lifting line method. Within the lifting line, the local loading distribution is calculated via the lift polar curve and the drag polar curve for each section of the wing (see Section 2.2.3). With STAR-CCM+®, it is necessary to cut slices of finite thickness in the wing geometry, in order to calculate the forces and torques on this slice. These are then divided by the thickness of the slice to get the local loading distribution. The mesh was refined around the slice (0.00625c r ) to minimize the irregularities in the mesh due to the slices (see Figures 3.4 A parametric study was needed to evaluate the influence of the slice thickness on the results (see Figure 3.5). For this purpose, three slices have been cut on the wing, one near the root of the wing, one near the tip and one in the middle of the half-wing. Thicknesses between 0.25c r and 0.001c r were tested. The converged thickness was estimated at e = 0.025c r . This thickness was used for all following simulations. RANSE simulations also allow the computation of the out-of-plane component of the local force on the section, which is not possible with the lifting line method. Simulations show that this component counts in average for 1.5% to 2.5% of the total force on the section.

2D results on the NACA2412 profile

For validation purpose, it has been decided to work with the NACA2412 profile instead of a LEI kite profile. Indeed the NACA2412 can be more accurately and easily meshed. Furthermore in following simulations, the profile was kept constant along the span of the wing.

2D RANSE simulations were carried out on the NACA2412 section for 14 angles of incidence within the range of [-8°; 16°]. Numerical results were compared with experimental ones obtained in a wind tunnel [START_REF] Abbott | Theory of wing sections[END_REF]) at the same Reynolds number, Re = 3.1 * 10 6 . The parameters of the parametric polar curves used for the lifting line method were chosen to fit the RANSE results at best. All these results are presented Figure 3.6. The agreement between experimental and numerical results is excellent for the lift coefficient and satisfactory for drag and moment coefficients. One explanation for the differences encountered for drag and moment is that simulations were conducted in a fully turbulent regime, whereas experiments were performed on a smooth section in a low residual turbulence wind tunnel.

Estimation of the numerical accuracy of the RANSE simulations

An attempt was done to estimate a numerical accuracy of 2D and 3D RANSE results. It was assumed that results from 2D simulations were sufficient for that purpose and that they could be extrapolated to the 3D cases. In addition, it was assumed sufficient to examine a single angle of incidence of 2 °for the root kite section, to be representative. size ranged from 0.6c r to 0.02c r , leading to a number of cells from 6 * 10 3 to 120 * 10 3 , keeping constant the near wall mesh parameters and the size of the computational domain (N = 10). Third (Figures 3.7 and 3.8(c)), the number of layers of the near wall mesh ranged from 4 to 32, leading to mean y+ values from 156 to 0.5, keeping constant all the other mesh parameters, that means the size of the computational domain N = 10, and the base size 0.1c r . In this last case, two other turbulence model were also tested (k-ω SST and Spalart-Allmaras).

For the variations of the domain size and of the number of cells (Figures 3.7 and 3.8, (a) and (b)), based on the limited number of computed points, for each aerodynamic coefficient, it was estimated a mean converged value (denoted cvg), a standard deviation relative to this mean value (denoted rms), and then a 95% confidence interval.

For the variations of the near wall mesh and of the models (Figures 3.7 and 3.8 (c)), the mean value of the coefficients from the different models (also denoted cvg) was estimated first, respectively around y+ ∼ 1 and y+ ∼ 30, then their standard deviations relative to these means due to model variations (also denoted rms). These are plotted in black Figures 3.7 and 3.8, and it was found that in each case the confidence interval encloses the red filled point corresponding to the general numerical set up. As a second step, the relative difference between coarse and fine mesh mean results (denoted dif ) was calculated. Estimated confidence intervals for the coefficients of lift, of drag and of moment, are resumed in the Table 3.1. Summing the first two and the last two columns of the Table 3.1, it was finally estimated that numerical results provided by the coarse meshed numerical set up of section were unconfined, mesh converged, and model independent, with the following relative accuracies: ±3.7% for the lift, ±7.7% for the drag and ±7.1% for the moment.

3D numerical results

For the validation of the lifting line method, 3D RANSE simulations have been carried out for some semi-circular wing geometries, from an un-twisted and un-swept wing with a linear chord law to a twisted wing with non-linear chord law and sweep law, the latter being close to a true kite geometry. The three different geometries are summarized in Table 3 These wings were set up at different angles of incidence and of sideslip. For the global aerodynamic results, the drag is computed as the component of the aerodynamic resultant in the direction of the relative wind, while the lift is computed as the component of the aerodynamic resultant in the normal plane to the relative wind. The aerodynamic moment is computed at the quarter chord point in the symmetry plane of the wing, along #» y b for the wing purely in incidence. For the wing in sideslip, the coefficients are computed from the norm of the aerodynamic forces. Global aerodynamic coefficients variations with respect to the angle of incidence or the angle of sideslip will be presented for each geometries. The error bars on these graphs are the estimated confidence interval at 95% for the STAR-CCM+® results (see Section 3.3) and for the mesh dependency of the lifting line method as described Section 2.3.2. Local results presented in following sections are the norm of the local force per unit length distribution along the span.

As an indication, the computation time for the non-linear lifting line method takes around a second on a desktop computer (Intel E3-1220, 4 cores, 3.10GHz, 4Go RAM) while each RANSE simulations lasted about 20 minutes on a cluster (Intel E5-2670, 8 cores, 2.6GHz, 128Go RAM).

Un-twisted and un-swept semi-circular wing

First, simulations have been carried out on a semi-circular with a 1.5m radius, untwisted and un-swept. The chord law varies linearly along the span from 1.0m at the root to 0.5m at the tips for an aspect ratio of 6.0. The results for this wing purely in incidence are presented in Figure 3.10. The difference between RANSE results and the lifting line method are around 5% for the lift coefficient. For the drag coefficient, there is a 5% gap at low angle of incidence up to 20% for the higher angles of incidence. The difference between the two approaches for the moment coefficient is of 12% in average over the investigated range.

Simulations were also carried out for this wing at 2°of incidence with different sideslip angles (see Figure 3.11). The difference is approximately of 10% for the lift coefficient, and of a few percent for the moment coefficient and the drag coefficient with a sideslip angle in the range of [0°; 15°]. At 15°of sideslip angle, results start to differ more significantly for the drag coefficient. As it can be seen, the loading and the moment on the wing increases considerably with the sideslip angle. For the use of the lifting line in a Fluid-Structure interaction, the local aerodynamic load per unit length has to be validated. In Figure 3.12, nondimensionalized magnitude of local forces per unit length are presented. Four simulations have been carried out, three at 2°of incidence with 0°, 7.5°, and 15°of sideslip (Figure 3.12 (a), (c) and (d) respectively) and another one at 12°of incidence only (Figure 3.12 (b)). Lifting line results follow the same trend as RANSE simulations and even in sideslip, minimum and maximum loading points are well predicted. The gap between the two methods is maximum for the maximum loading spot, and always lower than 15%.

Un-twisted and linearly swept semi-circular wing

To further validate the lifting line model, simulations have been carried out on a linearly swept wing with an aspect ratio of 7.5. The geometry is semi-circular with a 1.5m radius, un-twisted, with a constant sweep law of 30°along the span of the wing. This creates a very sharp angle at the center of the wing, thus a conical mesh refinement was added to enhance the simulation of the flow around the symmetry plane (see Figure 3.13).

Figure 3.14 shows the lift, drag and moment coefficient for the kite purely in incidence. The relative difference between the RANSE results and the lifting line model are around a few percent for the drag coefficient and 9% for the moment coefficient.

For the lift coefficient, the relative gap is around 7% and the difference between the two slopes of the lift coefficients in the linear range is of 10%. Local aerodynamic forces have also been computed for the wing at -4°, 0°, 10°and 15°of incidence, see Figure 3.15. The curves show a good consistency at the sides of the kite in contrast with the center of the wing, where local forces can differ of almost 60%. These results have to be put in perspective with the fact that the sweep angle is high (30°), therefore the angle in the middle of the wing is really sharp. This local singularity of the variation of the geometry induces locally a strong increase of the vorticity of the wake. This is well captured by the lifting line model which is intrinsically non diffusive. On the contrary, for the RANSE simulations, the mesh is coarse (y + = 35), even with the addition of the conical mesh refinement in the symmetry plane. The RANSE simulations are too diffusive to resolve the local wake vortex, and thus the associated loading variations. The local loading differences of sharp angle effects on the flow.

One simulation on this wing (case 2) with a sideslip angle has been carried out. The kite is at 5°of incidence with a sideslip angle of 7.5%. For the global coefficients, the gap between RANSE and lifting line results is of 19% for the lift coefficient and less than 6% for the drag coefficient. The local loading has also been computed on this swept wing and is shown in Figure 3.16 alongside the case at 5°of incidence and 0°of sideslip. Even with the difference in the middle of the wing, the position of minimum and maximum loading are still correctly estimated for the case in sideslip.

These results are obviously less satisfying than the first test case (case 1, un-swept and un-twisted wing) even if global results were still well predicted. Eventually a last investigation has been performed on a more realistic case, closer to a kite geometry, with a non-linear variation of the sweep law and without discontinuity in the symmetry plane.

-

Linearly twisted and non-linearly swept wing

The wing has a circular geometry with a 1.0m radius for an aspect ratio of 5.2, and a non-linear swept law (maximum swept at the tips of the wing : 0.5m). The chord law is also non-linear along the span and varies between 1.0m at the root and 0.1m at the tips. This wing is twisted linearly from 0°at the root to 5°at the tips (see Figure 3.17). This last wing is close to the geometry of a commercial leisure kite.

As previously, simulations were carried out for different angles of incidence ([-4°; 15°]) and for two angles of sideslip at 5°of incidence (see Figure 3.18). Purely in incidence, the lift difference stays below 5% and and moment coefficient difference is around 20%. The drag coefficient difference at low angles is only of a few percent, up to 30% at 15°of incidence. With an angle of sideslip, the difference is of 8% for the lift, from 2 to 15% for the drag and roughly 10% for the moment coefficient. Once more, the total loading of the wing increases with the sideslip angle and the variation in moment is a far more significant function of the angle of sideslip than the angle of incidence.

Local aerodynamic results are presented in following figures, in Figure 3.19, the wing was at -4°, 0°, 10°and 15°strictly in incidence. Once again, lifting line results follow the same trend as STAR-CCM+ ® . The gap between the two methods at -4°and 0°of incidence is larger than at higher angles of incidence but this has to be put in perspective since the global force is also lower than for the high angles of incidence, therefore the estimation error is eventually less significant than the one at higher angles of incidence. Furthermore, these low angles of incidence are not the usual flying angles of kites (see [START_REF] Borobia | Flight-Path Reconstruction and Flight Test of Four-Line Power Kites[END_REF]; Williams et al. (2008b)).

Finally, Figure 3.20 shows the local results on a wing at 5°of incidence and 0°, 7.5°and 15°of sideslip (respectively Figure 3.20

Conclusion

As shown previously, the global results of the 3D non linear lifting line method are really satisfying in comparison with the RANSE simulations, despite a very different approximation level, and over a large range of angle of incidence and of sideslip angle typically included between 0°and 15°. For the tested geometries, and especially for the one close to a LEI kite wing, lift coefficient estimations have less than 5% difference. The gap for the drag coefficient is a bit more significant for angles of incidence greater than 10°but this represents only a small portion of aerodynamic forces prediction on the wing. In sideslip, the lifting line is able to predict correctly the evolution of aerodynamic coefficients. Furthermore, the gap between the two methods is roughly maintained constant between the kite purely in incidence and the kite with a sideslip angle, which means that the error in incidence is kept but not increased with a sideslip angle. Local results are also satisfying, the lifting line is able to predict correctly the position of the minimum and maximum loading along the span, even in sideslip.

4 Fluid model application

Introduction

The modeling of a kite in flight is a challenging research topic as a kite develops its maximum power during a dynamic flight thanks to its own velocity. Most of the time the problem is studied from an aerodynamic point of view, with a rigid body kite, for simplification purposes. Leloup [START_REF] Leloup | Model approach and numerical tool development for kite performance assesment and mechanical design; application to vessels auxiliary propulsion[END_REF]) proposed an iterative procedure based on the zero mass model to calculate the effective value of the lift-to-drag ratio at each point of a given trajectory. Nevertheless, the kite is always assumed to be in pure translation against the apparent wind, and its turn rate is neglected, while it certainly induces an asymmetrical loading and may change significantly the equilibrium.

Williams (Williams et al. (2008b)) uses on the contrary a point mass model of the kite in order to determine the optimal trajectory to maximize power generation while keeping low the efforts needed to control the kite. Here again, the tethers are supposed straight but mass and drag effects are taken into account, assuming that the aerodynamic forces tangential to the line are negligible. The computation of the drag forces are discretized to consider the variation of movement speed, and thus wind velocity along the tethers. The aerodynamic forces on the kite are this time a linear function of the angle of attack and roll angle of the kite. However, as the kite is modeled as a point, its turn rate is yet again not considered.

In this chapter, a new procedure is proposed to calculate the effective lift-to-drag ratio of a kite along a prescribed flight path, taking into account the turn rate influence by means of the 3D Non-Linear Lifting Line model (LP3DNL). In a first step, the investigated trajectory is a circular trajectory around the full power point. Due to symmetry consideration, the flight of the kite is steady in this case and its turn rate can easily be calculated.

Kite in circular flight

In this model, the kite is reduced to the point K, and the tethers are reduced to only one rectilinear equivalent tether, whose direction and constant length are # » z k0 and L T , such that

# » OK = -L T # »
z k0 , O being the anchor point (see Figure 4.1). The weights of the kite and of the tethers are neglected. The trajectory is circular around the full power point, (O, # »

x RW ); the wind speed does not depend on altitude. The wind V RW # »

# » V RW # » z RW # » x RW # » y RW # » Va # » x V K # » x k0 # » y k0 # » z k0 R O K #» x b #» y b #» z b
x RW , the length L T of the equivalent tether, the radius R of the trajectory of the point K, and the kite geometry being prescribed, the goal in this approach is to find the equilibrium position of the kite, still assuming zero mass and the relative wind

# » V a = V a # »
x a parallel to the symmetry plane (K, #»

x b , #» z b ). The vector # » x k0 is defined as the direction vector of the projection of # »

x a or #» x b on the tangent plane in K at the flight sphere. The basis of the tangent plane is completed with # »

y k0 = # » z k0 × # »
x k0 . Three pseudo-degrees of freedom are introduced to set the kite relatively to the tangent plane in K at the flight sphere: a geometrical angle of incidence α g around (K, # » y k0 ) from #»

x b to # » x k0 , a roll angle γ around (K, # » x k0 ) from -# » y k0 to the projection of #» y b , and a yaw angle β around the tether axis (K, # » z k0 ) from # » x k0 to -# » x V K . α g is given, is kept fix, and is supposed to result from an adjustment of the real front/back tethers. γ and β are unknown variables. γ is supposed to result from an adjustment of the real left/right tethers, while β is supposed to result from a self-positioning of the kite.

Since this flight case is steady, the choice of point K along the studied circular trajectory doesn't matter, the flight direction # »

x V K is known, and the iterative algorithm starts with an estimated lift-to-drag ratio (0) and a roll angle γ (0) = 0. Denoting the counter of the recursive steps p, the kite velocity V (p) K can be calculated via the zeromass model formula with Equation 4.1 proposed by [START_REF] Leloup | Model approach and numerical tool development for kite performance assesment and mechanical design; application to vessels auxiliary propulsion[END_REF]), then the apparent wind by

V (p) a # » x a (p) = V RW # » x RW -V (p) K # » x V K .
The rotation rate of the kite is also calculated with Equation 4.3.

V (p) K = V RW   # » x RW • # » x V K + # » x RW • # » x V K 2 + # » x RW • # » z k0 sin( (p) ) 2 -1   (4.1) V (p) K = V RW # » x RW • # » z k0 sin( (p) ) 2 -1 = V RW 1 -R 2 L 2 T sin 2 ( (p) ) -1 (4.2) #» Ω (p) = V (p) K R # » x RW (4.3)
By the projection of p) on the tangent plane, # » x k0

# » V a (
(p) and β (p) are obtained. Using the 3D Non-Linear Lifting Line model for the kite in the wind # » V RW at α g , γ (p) , and β (p) , its kinematic defined by

# » V K (p) and #» Ω (p) , one obtain the aerodynamic resultant # » F a (p)
and the updated lift-to-drag ratio (p+1) . This force is naturally misaligned with respect to the equivalent tether direction # » z k0 , while it should not be the case according to the zero-mass hypothesis. To partially satisfy this equilibrium condition, the roll angle is firstly updated to γ (p+1) in order to get

# » F a (p) in the (K, # » x k0 (p) , # » z k0 )
plane. This has the effect of getting the apparent wind velocity out of the symmetry plane of the kite

(K, #» x b , #» z b ).
To recover this condition, the yaw angle is secondly updated to β (p+1) . The whole process is then iterated until convergence of , γ, β, according to some small numerical tolerance of 10 -8 .

It is found that this algorithm converges rapidly, and requires between approximately 5 iterations in the case of trajectories of big relative radius R/L T (leading to small roll angles), and 20 iterations for small relative radius cases (leading to high roll angles). When convergence is reached, the expected equilibrium is obtained and characterized by: the effective lift-to-drag ratio, the remaining attitudes angles γ and β, the aerodynamic resultant being effectively aligned with the equivalent tether and the apparent wind effectively parallel to the symmetry plane.

Results

Numerous calculations have been carried out in order to evaluate the variations of the lift-to-drag ratio , of the aerodynamic force # » F a , and of the angles of yaw β and of roll γ, with respect to the radius of the trajectory R, the tether length L T or the wind velocity V RW . A fixed angle of incidence α g of 3°have been chosen for the kite. Figure 4.2 shows the influence of the turn rate on the local loading of the kite. The kite follows a trajectory of radius 10m, with 50m of tether length, in a wind velocity of 5m/s. At convergence, the kite has a velocity of 51m/s for a lift-to-drag ratio of 10.41, with the lift estimated from the aerodynamic resultant taken orthogonal to the apparent wind direction. The black curve represents the relative local loading on the kite as obtained by the iterative equilibrium procedure described section 4.2, while the gray one shows the local loading of the kite with the same apparent wind but assumed to be purely in translation, the turn rate being neglected. The asymmetry of the loading caused by the induced velocities by the turn rate is clearly highlighted. 4.3 (a) shows the evolution of the lift-to-drag ratio depending on the relative radius, for different lengths of the tether between 25m and 200m. The lift-to-drag ratio is quasi-constant for longer tether lengths, in contrast it varies substantially for the smaller tether lengths. The same kind of variations can be found for the coefficients of lift and of drag (Figure 4.3 (b) and (c)). However in this particular case, a diminution of the lift-to-drag ratio corresponds to an increase of the lift coefficient and a faster increase of the drag coefficient. tether length, the kite velocity is almost the same, while it is more disparate for small relative radius (Figure 4.4 (b)). Nevertheless, despite the velocity decrease, the turn rate stays high due to the small radius of the trajectory.

V RW 5m/s C L 0.52 L T 50m C D 0.05 γ 14.6°V K 51m/s C L = L 0.5ρV 2 a S C D = D 0.5ρV 2 a S

Discussion

A first analysis have been made to study the influence of the real wind velocity V RW by keeping constant all other parameters. By dimensional analysis, the lift-to-drag ratio, the roll and yaw angle do not vary with respect to V RW while the kite velocity V K , its turn rate Ω and the apparent wind V a depend linearly on V RW . This derives logically from the linearity of Equation 4.1 and Equation 4.3 (see Section 4.2).

As it can be seen in Figure 4.3 (a), the lift-to-drag ratio tends to a constant value when the radius of the trajectory increases, the tether length being fixed. The relative radius being fixed, the limit values also reach a constant, when the tether length is increased. This constant value is the lift-to-drag ratio which is obtained neglecting the turn rate influence. Indeed, by increasing the radius or the relative radius of the trajectory one comes closer to cases of pure translation motion.

The lift-to-drag ratio decreases for small radii of the trajectory, while the lift and drag increase (Figure 4.3), even with a drop of the kite velocity (Figure 4.4 (b)). This augmentation of the aerodynamic forces is due to the high turn rate of the kite for the small radii of trajectory, which induces high local effective velocities. For the higher radius over tether length ratio, the lift-to-drag ratio is more important for the small tether length. For these cases, the drag force is quasi-identical (between 0.041 and 0.0460) while the lift force differs more (between 0.463 and 0.486).

The kite velocity is quite significant. It should be noticed that the equivalent tether drag is not taken into account in the iterative equilibrium solution presented in section 4.2. The following values are given as an example for the case recalled in Figure 4.2. The tethers drag has been calculated as Argatov et al. [START_REF] Argatov | Apparent wind load effects on the tether of a kite power generator[END_REF]) and Bigi [START_REF] Bigi | Investigation of the dynamic motions and operability of a ship towed by kite[END_REF]) by using the formulas from Hoerner [START_REF] Hoerner | Fluid-dynamic drag[END_REF]) with four lines, each with a diameter of δ = 2mm. The new lift-to-drag ratio is 9.13 instead of 10.41, which reduces the kite velocity to 44.7m/s instead of 51.0m/s, after convergence. Furthermore, the results are coherent with Dadd [START_REF] Dadd | Determination of kite forces using three-dimensionnal flight trajectories for ship propulsion[END_REF]) and Loyd [START_REF] Loyd | Crosswind kite power (for large-scale wind power production)[END_REF]), who stated that the velocity of the kite near the downwind position is close to the wind velocity times the lift-to-drag ratio of the kite.

The aerodynamic moment is also not taken into account. For the case recalled in Figure 4.2, the aerodynamic moment at point K is non null, but represents only a small part of the aerodynamic forces on the kite. The pseudo-roll component along # » x k0 is found to be the major one, the pseudo-pitch component along # » y k0 is the medium one (5 times smaller), and the pseudo-yaw component along # » z k0 is the minor one (15 times smaller). This last result for the yawing moment is consistent with the assumption made in section 4.2, which states that the yaw angle only results from a self-positioning of the kite.

Structure models presentation 5.1 Introduction

A leading edge inflatable (LEI) kite is made of an inflated tube as leading edge, a canopy and several battens to shape the canopy (see Figure 5.1). As a flexible structure, the modeling of the kite is a complex structural problem where optimal design still remains a crucial need. 2010)) developed a multi-body model to simulate the dynamics of the kite. The inflatable beams are modeled as a chain of rigid elements, linked by spherical joints which allow rotation along the three axis. The behavior of the inflatable beams is modeled by a torque vector on each joints. The torque vectors values come from experimental data and depend in particular of the internal pressure of the tube.

The kite canopy is modeled as a chordwise wire composed of discrete elements with an infinite stiffness. Those elements are linked together with hook joints which allow all rotations except twisting. Eventually the shear modulus of the fabric is modeled by springs which connect the nodes of the canopy but their stiffness do not depend on the aerodynamic pressure. The model is able to show global and local deformations, in particular a skewing deformation along the global X axis.

An even more complex model is described by [START_REF] Bosch | Finite element analysis of a kite for power generation[END_REF]; [START_REF] Bosch | Dynamic Nonlinear Aeroelastic Model of a Kite for Power Generation[END_REF]), with a finite element discretization of the wing. The canopy is modeled by nonlinear triangular shell elements, the beams by a linear Bernoulli beam element. In this model, the canopy is attached to the midpoint of the beam instead of its top. Bosch used a very coarse mesh and do not model small-scale phenomena such as wrinkling in the fabric because of computation cost.

For the mechanical properties of the fabric, tensile tests at 0°(warp direction), 45°and 90°(weft direction) allowed the identification of the orthotropic membrane stiffness matrix (in N/m) without measuring the thickness of the fabric, which is difficult to do. The tests have been carried out by Porcher Industries (a company member of the consortium) and by the laboratory 3SR in Grenoble, France. The longitudinal modulus E mc and E mt , the Poisson ratio ν ct , the shear modulus G mct , the longitudinal linearity limits X mc and Y mt and in shear S mct have been determined, where subscript m indicates the membrane aspect, c for the warp direction and t for the weft direction.

Inflatable beam modeling

The modeling of the inflatable beams, leading edge or battens, is quite significant for the global stiffness of the kite. Indeed, inflatable tubes serve both as a buoyancy for the kite and as stiffeners even without aerodynamic loading. As a structure exposed to significant stresses in various directions, the kite can be subjected to local buckling phenomena. A large bending due to a heavy aerodynamic loading for example, can trigger a sudden buckling of the leading edge beam initiated by local wrinkling of the fabric (see Figure 5.2). This buckling phenomenon on inflatable beams has been increasingly studied as the use of inflatable tubes for structural applications have been growing. Veldman [START_REF] Veldman | Design and Analysis Methodologies for inflated Beams[END_REF]; [START_REF] Veldman | Bending of anisotropic inflated cylindrical beams[END_REF]; [START_REF] Veldman | Wrinkling prediction of cylindrical and conical inflated cantilever beams under torsion and bending[END_REF]) investigated the wrinkling of inflatable beams under torsion and bending loads with experimental and finite element analysis. Veldman shows that a thin shell element can be preferred to a membrane element for the modeling of the fabric of the tube depending on the internal pressure of the inflatable beam.

Le [START_REF] Le Van | Bending and buckling of inflatable beam: some new theoretical results[END_REF]) developed a model with linearized equations from the Timoshenko beam theory for the study of the bending and buckling cases of an isotropic inflated beam. The analytical solution was then compared to finite element analysis, the inflatable tube being modeled by membrane elements. The paper [START_REF] Behrel | Investigation of kites for auxiliary ship propulsion: experiment setup, trials, data analysis and kite specs novel identification approach[END_REF] shows the influence of internal pressure on mechanical properties of the beam and highlights the better results of the Timoshenko beam theory compared to the Euler-Bernoulli beam theory for the modeling of inflatable beam.

Apedo [START_REF] Apedo | Nonlinear finite element analysis of inflatable beams made from orthotropic woven fabric[END_REF]) presents two finite element models (a linear and a nonlinear) of an inflatable orthotropic woven fabric beam. These two models are compared with a finely meshed thin-shell model and show a good agreement. Consequently beam models can be used for the modeling of inflatable beam if there is no wrinkling.

However, since the aim of this thesis is not to develop a new beam model, the threedimensional linear Timoshenko beam of Abaqus™ (B31 element) will be used. Due to this modeling, the internal pressure of the inflatable beams is not taken into account.

Nevertheless, the internal pressure has been used in [START_REF] Maison | De la problèmatique du dimensionnement de grands kites[END_REF]) in post processing, in order to determine the stress state in the tube wall.

The properties of the inflatable tubes are described by following notations: R B denotes the beam radius, E mc and G mct denote the mechanical properties of the fabric.

Following Expressions 5.2 to 5.9 were introduced in order to determine the different elementary section stiffnesses of the tube (see Equation 5.1):

          
Elongation stiffness:

E B A 0 = 2πR B E mc Bending stiffness: E B I = πR 3 B E mc
Transverse shear stiffness of the section:

G B (0.53A 0 ) = k Torsion stiffness: G B J = 2πR 3 B G mct (5.1)
The section area of the beam is set to:

Modeling of the canopy

The second part of the modeling of a LEI kite focuses on the canopy. To model such a piece of fabric, one can draw from sail modeling (see Augier et al. and Trimarchi et al. Augier (2012); [START_REF] Augier | Dynamic behaviour of a flexible yacht sail plan[END_REF][START_REF] Trimarchi | The use of shell elements to capture sail wrinkles, and their influence on aerodynamic loads[END_REF]). The canopy has the mechanical behavior of a membrane element, with a very small bending and transverse shear stiffness. Therefore modeling the canopy with a membrane element could be a logical choice. But a membrane has no bending stiffness and no compressive strength [START_REF] Lu | Finite element analysis of membrane wrinkling[END_REF]), which means that buckling or wrinkling can happen under a very low compressive load level. With a bending stiffness equal to zero, it is really difficult to predict the wrinkling or buckling phenomenon with a membrane element. However, it is possible to add a wrinkling model to the formulation of the membrane element (see Jarasjarungkiat et al., Schoop et al. and Barsotti et al. Jarasjarungkiat et al. (2008); [START_REF] Schoop | Wrinkling of nonlinear membranes[END_REF]; [START_REF] Barsotti | Static response of elastic inflated wrinkled membranes[END_REF]).

An other option for the modeling of this membrane behavior is the use of a thin-shell element. The shell element, with 5 degrees of freedom per node when the membrane element has only 3 per node, is generally more complex and computationally expensive (see Onate et al. [START_REF] Oñate | Membrane Structures Formed by Low pressure Inflatable Tubes[END_REF]). However, the computation is also more stable (see Trimarchi et al. and Jetteur at al. Trimarchi et al. (2010); [START_REF] Jetteur | Advanced Capabilities for the Simulation of Membrane and Inflatable Space Structures Using SAMCEF[END_REF]) since the missing of of bending and transverse shear stiffness creates some numerical instabilities. Trimarchi, Jetteur and Onate show the interest of thin-shell element to model a membrane structure, from an inflated tube to a spinnaker sail.

Because numerical instabilities and numerical convergence were a real problem when using membrane element, notably with the first Kite as Beam model presented in Solminihac [START_REF] Solminihac | Kite as a Beam -An analytical 3D Kite Tether Model[END_REF]), the kite canopy is therefore modeled by the thin-shell element S4R5 of Abaqus™ . The Discrete Kirchhoff constraint is imposed on this shell element, which means that the Kirchhoff constraint is satisfied at discrete points on the shell surface ( see SIMULIA manual SIMULIA ( 2014)). To enforce this constraint, Abaqus uses the transverse shear stiffness, defined as:

Kts αβ = f p K ts αβ (5.10)
where Kts αβ is the component of the section shear stiffness (α, β refer to the local direction 1, 2 associated with the shell section), f p is a dimensionless factor used to prevent the shear stiffness from becoming too large in thin shells and K ts αβ is the actual shear stiffness of the section. f p is defined as:

f p = 1 1 + 2.5 × 10 -5 A t 2 (5.11)
where A is the area of the element and t is the thickness of the shell.

The goal of the following paragraph is to estimate a realistic value of the transverse shear stiffness of the shell element. The material is assumed to have an isotropic be-havior. The bending and transverse shear stiffness of the shell element is assumed to be equivalent to those of a Timoshenko beam of length L, height t and width b.

For the canopy, the shear stiffness of a beam element of length L is:

12EI L 3 1 + 12EI GA 0 L 2 (5.12)
As the transverse shear stiffness should be negligible compared to the bending stiffness:

12EI GA 0 L 2 1
(5.13)

The shear stiffness becomes for a rectangular section of height t and width b:

12EI L 3 = Eb t L 3 (5.14)
Assuming a square mesh, b = L, with t the thickness of the shell and A the area of the element:

Eb t L 3 = Et t 2 A
(5.15)

According to Equations 5.10, 5.11 and 5.15, the actual shear stiffness has to be set to:

K ts 11 = K ts 22 = Et t 2 A + 2.5 × 10 -5
(5.16)

Subsequently, the ratio between the membrane stiffness and the bending stiffness of the shell elements are also carefully chosen in order to obtain the same mechanical behavior between membrane and shell element and also to avoid numerical instabilities due to wrinkling. The numerical tests presented in Section 5.3.4 show that 10 -8 is an appropriate value for this ratio.

The bending stiffness D and the membrane stiffness C for an isotropic material behavior are defined as:

D = Et 3 12(1 -ν 2 ) and C = Et (1 -ν 2 )
(5.17) Which leads to following equation, for a given surface S of the modeled object, either a kite or an elementary cell (see Sections 5.2 or 5.3)

D CS = t 2 12S = 10 -8
(5.18)

The thickness t of the shell element is then defined in relation to the surface of the modeled object; such as

t = √ 12S10 -8
(5.19)

Complete Finite Element Model

This model is an evolution of the Finite Element Model presented by [START_REF] Leloup | Model approach and numerical tool development for kite performance assesment and mechanical design; application to vessels auxiliary propulsion[END_REF]). This Section includes a lot of developments previously done by [START_REF] Maison | De la problèmatique du dimensionnement de grands kites[END_REF]) who worked on the beyond the sea project at the same time of this thesis. In the rest of this document, we will refer to this model called the IDK model, acronym for the french abbreviation "Intégration fluide-structure Dimensionnement Kite".

As the IDK model is quite heavy, only the half of the kite is modeled. This reduces the simulation possibilities to symmetric cases but on the other hand also decreases the computation time. The chosen Abaqus™ solver is based on an implicit dynamic method with a quasi-static application (backward Euler time integration). This method is used for nonlinear quasi-static applications on structures of small stiffness and large displacement. It avoids convergence issues since high numerical damping is produced with this integration scheme. With this method, it is necessary to carry out independently a study of the eigen frequencies in order to determine the lowest eigen frequency since the total time of the dynamic simulation should be at least slightly above the equivalent period. In order to confirm the correct convergence of the simulation, it is necessary to check the global kinetic energy which should be negligible in comparison to the global strain energy. If not, the total time step should be increased.

The shape of the canopy is defined from the LEI kite profile determined by [START_REF] Leloup | Model approach and numerical tool development for kite performance assesment and mechanical design; application to vessels auxiliary propulsion[END_REF]) and presented in Figure 5.3. In the IDK model, it is assumed that the same dimensionless profile is used along the span of the wing. As the inflatable leading edge is modeled as a beam, the beam elements are put at the end of the canopy, on the black cross of the Figure 5.3.

The radius of the leading edge beams vary linearly from the middle of the kite to the tips as a function of the span of the wing. The inflatable battens are placed on the canopy mesh. Thus, they have the same shape than the canopy, that means the LEI kite profile. The radius of the inflatable beams varies also linearly between the leading edge and the trailing edge of the wing. For the loading of the model, a pressure can be applied on the shell element, which accounts for the canopy. Forces and moments due to pressure integration on the 3D leading edge can also be applied on the leading edge beam elements.

Kite as a Beam

The Kite as a Beam (KaB) described here is an evolution of the model presented by de Solminihac [START_REF] Solminihac | Kite as a Beam -An analytical 3D Kite Tether Model[END_REF][START_REF] Solminihac | Airborne Wind Energy[END_REF]). The core idea of the Kite as a Beam model is to approximate a LEI kite by an assembly of equivalent beam elements (see Figure 5.5). Only isotropic behavior of the canopy material is implemented in this version, but it can be upgraded to orthotropic formulation with some adjustments. For the isotropic modeling of the fabric, the mean value of Young modulus and Poisson ratio in warp and weft directions is used. The wing is partitioned into several elementary cells (see Section 5.3.1), each one modeled with a Timoshenko's beam of equivalent mechanical properties (see Section 5.3.2). Eventually, the equivalent beams are re-discretized in order to match with the lifting line mesh (see Section 5.3.3). The structural analysis is performed with the finite element solver Abaqus™.

Elementary cell

Inflatable batten Equivalent beam

The elementary cell concept

The concept of the Kite as a Beam is to model an LEI kite as an assembly of equivalent beams. The kite is discretized in several elementary cells, which consist of two inflatable battens, the leading edge tube in between and the corresponding canopy. In this model, the tips of the wing are considered as battens (as shown in Figure 5.5).

To simplify the computation of the mechanical properties, the elementary cell is approximated by a rectangle and its four corners belong to the same plane. The spanwise width L LE and the chordwise length L BA of the cell are extracted from the geometry of the kite (see Figures 5.6 and 5.7).

A L A R F R F L M R M L L LE #» x c #» y c Figure 5.6 -Elementary cell geometry # » M R = 0.75 # » A R + 0.25 # » F R # » M L = 0.75 # » A L + 0.25 # » F L
(5.20)

L LE = # » M R M L (5.21) L BA = 1 2 # » A R F R + # » A L F L (5.22)
The local coordinate system of the cell

( #» x c , #» y c , #» z c
) is calculated as follows:

#» x c = # » M L M R L LE (5.23) #» z c = # » A L F R × # » A R F L # » A L F R × # » A R F L (5.24) #» y c = #» z c × #» x c
(5.25)

Eventually, the geometry of the canopy in the elementary cell has to be determined. Indeed the geometrical stiffness of the canopy can not be neglected and its initial shape has an influence on the total stiffness of the cell (see Section 6.3). The aim of the technique is to be close to the geometry of the real kite but also to keep a straightforward geometry construction. To this end, the battens are straight and the canopy is curved following #» x c , as it can be seen in Figure 5.7. The curvature is the sum of a 2nd degree polynomial according to #»

x c and a sine according to #» y c (see Equation 5.27). The 2nd degree polynomial represents the distance between the real geometry of the kite and the plane of the cell and the sine represents the camber of the Leading Edge Inflatable kite section. The geometry is given by Equations 5.26 to 5.28 with x ∈ [0; L LE ] and y ∈ [0; L BA ]. L HT is defined as the distance between the middle of the segment M R M L and the middle point of the lifting line included between M R and M L . Knowing the geometry of the kite section profile, C MX represents its camber and R LA the relative position of the maximum camber (see Figure 5.7). The following values C MX = 0.065 and R LA = 0.3 allow the description of a canopy adequately close to the LEI kite section. Each point of the surface thus created can be defined by the vector position #» x according to the following expression:

#» x = x #» x c + y #» y c + z #» z c
(5.26)

with z = - 4L HT L 2 LE (x - L LE 2 ) 2 + L HT + CL BA sin y L BA γ π (5.27) where C = x 4C M X L LE 1 - x L LE and γ = log(1/2) log(R LA )
(5.28)

M L M R L HT L HT + C MX L BA R LA L BA L LE # » xc #» yc #» zc L BA

Shell element

Lifting line Beam element Truss element As it can be seen in Figure 5.7, the inflatable battens (leading edge and battens) are modeled by linear beam elements. Starting from the known beam radius R B , the fabric properties E mc and G mct , the section properties are estimated as described in Section 5.1.1.

Because in the complete wing model the elementary cells are connected, the following linear superposition assumption is made: the stiffness of the finite element beam representing a batten between two cells is only 50% of the stiffness of the full batten. For the battens at the tips of the wing, the stiffness of the finite element beam is the full stiffness of the batten.

For the modeling of the canopy, shell elements have been chosen (see Section 5.1.2). The Young modulus E CA , the density and the Poisson ratio are known. The thickness of the shell element e CA depends of the geometry of the cell (see Equation 5.19):

e CA = 12L BA L LE 10 -8
(5.29)

K 11 = E CA e CA e 2

CA

A elem + 2.5e -5

(5.30)

The last part of the modeling of the elementary cell is the leech line on the trailing edge of the wing. This is modeled by truss elements with Dyneema properties: the density ρ DY and the Young modulus E DY .

Both the leading edge and the inflatable battens are discretized with 200 beam elements and the canopy with about 4000 shell elements (see Section 5.3.4).

The equivalent beam properties

To compute the equivalent beam properties, the elementary cell is set under pressure and then subjected to various linear perturbations. Knowing the reaction forces to these linear perturbations allows the determination of the equivalent beam mechanical properties. Many options can be used to estimate the equivalent beam properties.

In particular, several elementary cases of loading and associated boundary conditions have been tested. The model presented in this subsection has been selected considering that the global opening (distancing between the two opposites tips) and twist of the kite are the most influent deformations that govern the aerodynamic loading.

The first computation step is a non-linear dynamic analysis for quasi-static application.

The four corners of the elementary cell are clamped and the cell is loaded with an estimated homogeneous pressure. Each cell can have its own pressure. Equation 5.31 gives a good estimation of the pressure, with ρ the air density, C L the lift coefficient of the wing for the respective apparent wind V a .

P 0 = 1 2 ρC L V 2 a (5.31)
A second non-linear static computation step is then conducted with the same conditions to check the validity of the solution obtained. A representative simulation result is shown in Figure 5.8.

Starting from this pressurized structure, six linear perturbation calculation cases were completed in order to evaluate the stiffnesses of the elementary cell with respect to x c , #» y c , #» z c ) in our model. For all cases, the beam B L (see Figure 5.8) is clamped while various boundary conditions are applied to each nodes of the whole beam B R . The results for the case (β) can be seen in Figure 5.9.

#» z c #» y c #» x c B L L e

Case

(α) (β) (γ) (δ) ( ) (φ) U      a 0 0           0 0 a           0 a 0      U 3 = 0 U 1 = 0 U 2 = 0 U R #» 0 #» 0 #» 0      0 0 ω           ω 0 0           0 ω 0     
Table 5.1 -Boundary conditions on the beam B R in displacements (U ) and rotations (U R) for the load cases (α)-(φ), components expressed in the frame

( - → x c , - → y c , - → z c )
Numerical results depend linearly on a and ω since a linear perturbation mode is used. Reaction forces and moments on the right beam B R are measured for each load case in the direction of the elementary displacement and rotation. F R(X) is the sum over all the nodes of the reaction forces RF on the right beam while M R(X) is the sum over all the nodes of the reaction moments RM .

In order to separate transverse shear and torsion, the distance of the median line of the equivalent beam from the leading edge is computed as:

D = M R(β) F R(β) with M R(β) = B R yRF 3(β) + RM 1(β) (5.32)
where y is the coordinate of the considered node, RF 3 is the node reaction force component along #» z c and RM 1 the node reaction moment along #» x c . The equivalent beam extremities P L and P R are given by the following equation, where M L and M R are the points of the lifting line:

P L = M L - L BA 4 -D #» y c and P R = M R - L BA 4 -D #» y c
(5.33)

The elongation stiffness is calculated as:

EA 0 = L LE a F R(α)
(5.34) and the torsional stiffness as:

GJ = L LE ω M R( ) (5.35)
The in-plane bending stiffness about #» z c is determined by:

EI 3 = L LE ω M R(δ)
(5.36) and the bending stiffness about #» y c by:

EI 2 = L LE ω M R(φ)
(5.37)

Abaqus™ computes the transverse shear stiffness as (see Equation 5.10):

K i1 = f i p K i1
(5.38) with

f i p = 1 1 + SCF l 2 A 0 12I j (5.39)
where i, j ∈ {2, 3} and j = i.

The Slenderness Compensation Factor SCF is set to 0.25 which is the default value of Abaqus™ and l is the length of the micro-beam element (a re-discretization of the equivalent beam to correspond to the mesh of the lifting line model, see Section 5.3.3 for more details).

Since the value of K 31 should be:

K 31 = L LE a F R(β)
(5.40)

The value of the transverse shear stiffness of a micro-beam element introduced in Abaqus™ is:

K 31 = 1 + l 2 EA 0 48EI 2 L LE a F R(β)
(5.41)

In the same way:

K 21 = 1 + l 2 EA 0 48EI 3 L LE a F R(γ)
(5.42)

Kite construction

The kite is now modeled as an assembly of several equivalent beams, whose mechanical properties are known. Since this model is supposed to be used in a Fluid-Structure interaction scheme, the properties of the equivalent beams are attached to the lifting line mesh (see Section 2.2.1). This means that the points M of the lifting line (see Section 2.2) define the tips of small finite element beams, called micro-beams (see Figure 5.10). The structural analysis of the kite is then conducted from this geometry.

Leading edge

Trailing edge Rigid body Micro-beam/Lifting line Equivalent beam Since the geometrical location of the finite element micro-beam lies on the lifting line, its local section direction # » n 1 is determined in the following manner. The point C is defined as the middle of the micro-beam element. The point P C is the orthogonal projection of C on the correspondent equivalent beam (see Figure 5.11).

P C = P L + ( # » P L C • #» x c ) #» x c (5.43) If #»
t stands for the unit vector along the micro-beam element axis, the unit vector # » n 1 is obtained from:

# » n 1 = # » CP C -( # » CP C • #» t ) #» t d mb (5.44) with d mb = # » CP c -( # » CP C • #» t ) #» t (5.45)
It is assumed that the location of the micro-beam element section centroid is expressed as d mb # » n 1 . The second local section direction of the micro-beam element # » n 2 is such that:

# » n 2 = #» t × # » n 1
(5.46)

The micro-beam element section properties are the same as defined in Section 5.3.2 assuming the micro-beam element frame (

#» t , # » n 1 , # » n 2 ) is matching the equivalent beam frame ( #» x c , #» y c , #» z c ). Figure 5
.11 is distorted for the explanation,

#» t is almost parallel to #»
x c . This assumption is discussed in Section 6.3.

For the loading of the model, forces and moments can be applied on the lifting line points, tips of the micro-beams.

Equivalent beam

Micro-beam

P L C # » xc P C P R #» t # » n 1 d mb

Numerical tests for the elementary cell

A mesh convergence study is needed to evaluate the influence of the elements number on the results. The mesh of the cell is controlled through N beam the number of beam elements. All the leading edge and the two inflatable battens are composed of N beam beam elements of equal size which generates a mesh with almost square element (as it can be seen in Figure 5.9). In Figure 5.12 and in Table 5.2, N beam ranges from 50 to 500 beam elements. Figure 5.12 represents the relative error (Equation 5.47) for the traction case (α) on the cell in the middle of the kite.

err = F R(α) (N beam = X) -F R(α) (N beam = 500) F R(α) (N beam = 500) * 100
(5.47)

Beams number

Relative error (%)

Figure 5.12 -Relative error on the (α) traction case with respect to the number of beam elements for the largest cell of the kite.

N beam 50 100 200 300 400 500 Total CPU time (s)

4.0 11.3 44.7 102.9 193.8 614.8

Table 5.2 -Computation time (total CPU time) of all the linear cases with respect to the beam numbers for the largest cell of the kite.

For N beam = 200, it is estimated that the results are converged enough (err = 0.16%). Furthermore, as it can be seen in Table 5.2, the computation time starts to increase significantly for N beam > 200.

As explained in Section 5.1.2, the ratio between the membrane stiffness time surface area C S and the bending stiffness D is crucial for the choice of the shell thickness. Figure 5.13 shows the relative error in the reaction forces of the (α) case for a ratio between 10 -5 and 10 -11 .

- Relative error (%) As the computation time increases to a not negligible extent with the ratio diminution since it also decreases the shell thickness, 10 -8 is chosen with a relative error of 5.94%.

Structure models comparison

The aim of this chapter is the comparison of two structure models: the Kite as a Beam (KaB) model and the more complex Finite Element model (IDK). The comparison cases considered are elementary displacements along the three global axes and a torsion of the wing tip, as the twist of the wing can be a significant deformation of the kite in flight condition. Comparison cases are detailed in Section 6.1, results are presented in Section 6.2 and discussed in Section 6.3

Description of the comparison cases

In order to correctly compare the two models, applied boundary conditions should be as consistent as possible, keeping in mind that the IDK model is a complete kite geometry and the KaB model is an 3D assembly of beams. As a working pressure has to be chosen to compute the properties of the equivalent beams in the KaB model (see Section 5.3.2), the IDK model should also be put under pressure. To compare the global stiffness given by the two models, they have to be applied for the same conditions. As the IDK model is deformed when put under pressure, the same deformations have to be applied to the beam assembly of the KaB model so that the two geometries are as close as possible. Both models are then successively submitted to the different elementary cases, using Abaqus™ .

Kite under pressure

As a first step, the IDK model is put under a uniform pressure (see Figure 6.1). For the computation, following boundary conditions are chosen. In the symmetry plane, the displacement along the X and Y axis are null as well as the angular velocity along the X and Z axis. At the tip, the displacement along the Z axis is canceled. A kinematic coupling is implemented at the plane of symmetry in order to suppress the rigid body translation of the kite along X, without over constraining the canopy. The point S, closest point of the quarter chord in the symmetry plane, is the reference point for this constraint (U X (S) = 0) which links each point of the symmetry plane.

At the end of the simulation, the positions of the leading and trailing edge are extracted. This allows the computation of the displacement # » D IDK of the equivalent lifting line geometry of the IDK model:

Trailing edge #» X #» Z #» Y Symmetry plane Wing tip U Y = 0 VR X = VR Z = 0 Kinematic coupling, U X = 0 at S U Z = 0
# » D i IDK = ( 3 4 # » A i def + 1 4 # » F i def ) -( 3 4 # » A i init + 1 4 # » F i init ) (6.1)
where i is the identification number of the meshed canopy, with A i the leading edge point, F i the trailing edge point, def and init stand for the distorted and the initial geometries respectively.

The next step is to apply the displacement calculated from Equation 6.1 to the mesh of the KaB model, that means the lifting line. The mesh of the KaB model is identified with the superscript k while the mesh of the IDK model is identified with the superscript i. As the two meshes do not correspond since the models do not have the same discretization, a correspondence between KaB points and IDK mesh nodes has to be found. Therefore, the point of the KaB trailing edge F k init is assumed to be between the points of the IDK trailing edge F i init and F i+1 init when d 1 < d and d 2 < d, as shown in Figure 6.2. The displacement to apply to the KaB geometry

# » D k
KaB is then computed as described in Equation 6.2.

# » D k KaB = ( d 2 d 1 + d 2 ) # » D i IDK + ( d 1 d 1 + d 2 ) # » D i+1 IDK (6.2)
It was chosen to work with the trailing edge to estimate the corresponding points of each geometry because it was easier to make them be coincident. The last step is eventually to re-estimate the orientation axis of the micro-beam function of the applied deformation. This is done by applying the computed displacement to the tips of the equivalent beam and then calculating the local frame of the micro-beam as described in Section 5.3.3.

F i+1 init F i init F k init d 1 d 2 KaB Mesh IDK Mesh d Figure 6
.2 -Details of the notations used for the distorted geometry equivalence between KaB and IDK

The four elementary comparison cases

Four test cases were chosen in order to evaluate the elementary deformation that can encounter a kite in flight. The first case is a displacement along the X axis of the wing tip. To evaluate the opening and closing of the kite, the second and third cases are respectively a displacement along the Y axis of the wing tip and along the Z axis of the kite root. The last case is the torsion of the wing tip (see Figure 6.3).

As the movement should be equivalent for the IDK and KaB models, a particular attention has been given to the boundary condition applied. For the IDK model, same boundary conditions as for KaB are propagated for the pressure step in addition to a velocity condition of the type:

V = a T (6.3)
with T the simulation time and a the chosen displacement for the prime three cases. Indeed a velocity condition was the best way to apply the same displacement on the set of points after the deformation of the pressure step.

For the torsion case, discrete displacement conditions are applied in order to impose a rotational motion. The rotation center C is the closest point from the quarter chord which is constrained along the Y axis. The displacement U Y imposed on a point P of the wing tip is calculated as follows:

KaB model wing tip for all the four cases. Table 6.1 resumes the boundary conditions for the IDK model.

Symmetry plane (S) Wing tip

X axis Boundary conditions for the KaB model are described in Table 6.2. With these conditions, it is assumed that the movement for each comparison cases are equivalent.

U X = U Y = 0 U Z = 0 VR X = VR Z = 0 V X = a/T Y axis U X = U Y = 0 U Z = 0 VR X = VR Z = 0 V Y = a/T Z axis U X = U Y = 0 VR X = VR Z = 0 U Z = 0 V Z = a/T Torsion U X = U Y = 0 U Z = 0 VR X = VR Z = 0 U Y (P ) = (P -C). #» X tan(ω)

Root (K) Tip (O)

X axis

U X = U Y = 0 U Z = 0 VR X = VR Z = 0 VR X = VR Y = VR Z = 0 U X = a Y axis U X = U Y = 0 U Z = 0 VR X = VR Z = 0 VR X = VR Y = VR Z = 0 U Y = a Z axis U X = U Y = 0 U Z = 0 VR X = VR Y = VR Z = 0 VR X = VR Y = 0 U Z = a Torsion U X = U Y = 0 U Y = U Z = 0 VR X = VR Z = 0 VR X = VR Y = 0 UR Z = ω

Comparison of the two models

For the analysis of the results provided by the two models, reaction forces at the wing tip along the X and Y axis are added for the X and Y cases respectively. For the Z case, the reaction forces at the symmetry plane along the Z axis are summed. For the torsion case, the reaction moment at the point O along the Z axis is an output of the KaB simulation. For the IDK model, the moment RM is computed from the reaction forces RF 2 along the Y axis and the distance projected on the X axis between the point P and the rotation center C (see Equation 6.5).

RM IDK = P (P -C). #» XRF 2 (6.5)
The difference between the two models is calculated as following with RF the reaction force result and RM the reaction moment:

error = 100 |RF IDK -RF KaB | RF IDK
for the translation cases (6.6)

error = 100 |RM IDK -RM 3 KaB | RM IDK
for the torsion case (6.7) Table 6.3 shows the average simulation time (on a classic computer) for the two models. As it can be seen in Table 6.3, the IDK computation time nearly doubles between a pressure of 100P a and 400P a while the increase is only of 25% for the KaB As it can be seen in Figure 6.5, the two models can differ significantly depending on the pressure. The IDK model varies quasi linearly with the pressure while it is not necessarily the case for the KaB model. For the Y and Z displacements, which characterize the closing and opening of the kite, the KaB reaction force is almost a constant function of the pressure. The gap between the two models are less than 50% except for the torsion case, where the KaB model overestimates largely the IDK results. But for the torsion case the majority of the deformation are at the tip of the wing, which is trapezoidal in the complex model and rectangular in the KaB model. This overestimation of the tip geometry may play a significant part in the gap between the two results.

Comparison with a cylindrical geometry

To evaluate the difference related to the geometric approximation of the KaB model, that is the elementary cell as a square, another set of simulations have been carried out with a basic geometry which can be seen in Figure 6.6. The kite is a half-cylinder, with a radius of 4m and a constant chord of 1m. There is no twist and no wingsweep. This wing geometry reduces the geometrical approximation when building the elementary cells to a minimum. The four corners of the cell already belong to a plan, each cell is already rectangular and without wingsweep the cell is perfectly oriented.

Simulation results can be seen in Figure 6.5 and the difference between the two models computed as described in Equations 6.6 and 6.7 can be seen in Table 6.4.

The stiffnesses given by the KaB model are still almost constant for the Y and Z cases, with a non-negligible difference compared to the stiffnesses obtained by the IDK model, except for the torsion and the X cases where the results are a bit closer. Results for the displacement along the X axis follow almost the same trend, with a gap lower than 25%. For the torsion case the gap is almost lower than 20% which is a Table 6.5 -Gap (in %) between the IDK and the KaB models for the cylindrical wing.

Discussion

Even when considering the level of approximation of the two models, there is still a significant gap between them. For the closing and the opening of the kite, i.e. the Y and Z displacement, the two models do not follow the same trend, even if they can coincide on some pressure values. For the X case, the results are more satisfying with a gap lower than 25% on both geometries. Eventually, the torsion is largely overestimated by the KaB model for the first geometry and this seems mostly due to the geometric approximations made when building the elementary cell. But the evolution of the KaB results for the cylindrical wing on the Y and Z axis shows that the difference is not just due to the geometric approximations.

For the computation of the properties of the equivalent beams, the boundary conditions are applied on the whole batten, meaning that the batten is assumed as a rigid body. The validity of this assumption has not been discussed before. Figure 6.8 shows the deformation of the battens in the IDK model under a pressure of 200P a and after an elementary displacement along the X axis. The batten deforms of about 3% at most from the geometry under pressure after displacement. Assuming the battens as rigid bodies does not seem to be the major explanation for the difference between the two models. A second explanation of the difference between KaB and IDK results is the description of the aerodynamic profile, that means the shape of the canopy. Indeed, the shape of the canopy in the KaB model is described analytically (see Section 5.3.1) as close as possible from the section of a LEI kite while the canopy of the IDK model is estimated directly from the aerodynamic profile of a LEI kite (see Section 5.2). As it can be seen in Figure 6.9, the difference may be significant depending on the influence of the canopy shape on the kite stiffness. .12 shows the equivalent beam properties for the cell with various positions of maximum camber (R LA ∈ [0.25; 0.5]) for the same maximum camber C MX = 0.065 under the tested range of pressure. The case R LA = 0.2 is not represented here since the sine discretization along #» y c creates a very deformed mesh around the leading edge which distorts the results of the simulation. The evolution of the properties are globally linear except for the in-plane bending stiffness about #» z c , EI 3 and the transverse shear stiffness about #» y c , K 21 . The curve slope of the stretching stiffness and the bending stiffness about #» y c (respectively EA 0 and EI 2 ) is quite small, which may explain the results of the KaB model for the comparison cases about the Y and Z axis, that means the small evolution of global stiffness function of the pressure.

There is almost no influence of the position of maximum camber on the torsional stiffness GJ, on the bending stiffnesses EI 2 and EI 3 , or on the stretching stiffness EA 0 and the transverse shear stiffness along #» z c , K 31 . However, the distance of the equivalent beam from the leading edge D increases with the rise of the coefficient R LA . The transverse shear stiffness along #» y c also vary with the increase of the coefficient R LA . Nevertheless the modification of the position of the maximum camber does not change drastically the curve of the mechanical properties, and will not change either the results for the comparison cases. Figure 6.13 shows the equivalent beam properties for the middle cell of the kite under a pressure from 100P a to 500P a, for a maximum camber C MX between 0.001 and 0.1, with the same position of maximum camber R LA = 0.3. As it can be seen in Figure 6.13, the properties are very similar for C MX values close to C MX = 0.065. Only C MX = 0.001 and C MX = 0.1 seem to differ significantly, in particular for the distance of the equivalent beam from the leading edge D. Figure 6.13 and 6.12 show that the shape of the canopy (if kept close from the LEI kite section) has only a small influence on the mechanical properties of the equivalent beam and thereafter on the difference between the IDK and the KaB models. Another major assumption of the KaB model is the equivalence between the equivalent beam frame ( #» x c , #» y c , #» z c ) and the micro-beam frame ( #» t , # » n 1 , # » n 2 ) as described in Section 5.3.3. Figure 6.14 shows the two frames in the middle of the kite, where the approximation is the largest. In the middle of the cell, #» t and #»

x c are equivalent but the angle difference between # » n 1 , #» y c and # » n 2 , #» z c is about 18°. On the contrary, around the battens, the largest difference is between #» t and #»

x c , and also about 18°for this LEI kite geometry. This hypothesis is inherent to the KaB model and the choice to match the structure and the lifting line mesh, and also to the description of the beam element in Abaqus™ .

Lifting line points Leading edge

Trailing edge Inflatable battens Micro-beam Equivalent beam Micro-beam axis Equivalent beam axis 7 Fluid-structure interaction

# » xc #» yc #» zc # » n 1 # » n 2

Introduction

Fluid-Structure interactions (FSI) is a growing field of research with the increase of available computational power. FSI are used in various domains, from civil engineering to biomedical (De Nayer ( 2008)), for compressible or incompressible fluids. This chapter will focus on FSI in the maritime industry, on interactions between incompressible air and membrane fabric such as sails. Indeed, sails problems are very similar to kite ones. The sail is a membrane structure with or without battens subjected to aerodynamic loadings and large displacements. As the performances of the sails depend on their shapes, FSI have been used since several years in sailing in order to predict the performances.

Charvet [START_REF] Charvet | Simulation numérique de l'écoulement autour des voiles de bateau[END_REF]) in 1991 uses a 3D particle method to determine the effects of the interaction between the jib and the main-sail and couples this method with a shell finite element model to take into account the deformations of the sail in a stationary case. Hauville [START_REF] Hauville | Dynamic control of a sail by fluid/structure intercation method[END_REF]; [START_REF] Hauville | Equilibre dynamique d'une structure idéalement flexible dans un écoulement : application à la déformation des voiles[END_REF]) uses a wire grid to represent the sail, the warp and the weft of the fabric. This model is coupled with a particle method which allows the computation of the strain on the sail. As the structure mesh is triangular while the fluid mesh is quadrangular, the pressure on each node of the structure mesh is estimated with a bi-cubic interpolation. The iterations between the two models continue until an equilibrium is found. The dynamic effects can also be estimated with an iterative method, function of time.

Malpede [START_REF] Malpede | A fully integrated method for optimising fibermembrane sails[END_REF]) also worked on sails and developed a method to find the optimal fiber layout distribution in fiber membrane sails. As geometric nonlinearities can occur due to large displacement, the structure model uses a nonlinear finite element method even if the material properties are kept linear. An adaptation of a vortex lattice aerodynamic model is used for the aerodynamic analysis. Here again, the method iterates between structure and aerodynamic models and a bi-cubic interpolation method is used to apply the pressure load on the finite element nodes.

In the last years, as the use of kites for energy production has developed, Fluid-Structure Interactions on kites became a research subject, which is still under the highlights. Breukels [START_REF] Breukels | An engineering methodology for kite design[END_REF]; [START_REF] Breukels | Airborne Wind Energy[END_REF]) uses a multi-body system to model the kite with canopy and battens, and the tethers (see Section 5.1). Breukels carried out CFD simulation on a 2D airfoil to compute its lift, drag and moment coefficients as a function of the angle of incidence, thickness and camber of the airfoil.

The pressure distribution can also be determined with the CFD simulation, which gives the forces distribution on the nodes of the canopy. The section angle of incidence are calculated with a vortex lattice method for the kite wing. Breukels compared a model simulation with experimental data, tension in the four lines of various kites during a loop. The velocity of the kite and the wind speed were also measured. Measured forces on the steering line are used as control input for the simulation, while the forces in the power lines and the velocity of the kite are compared to the experiment results.

Bosch [START_REF] Bosch | Finite element analysis of a kite for power generation[END_REF]; [START_REF] Bosch | Dynamic Nonlinear Aeroelastic Model of a Kite for Power Generation[END_REF]) is an other example of FSI on kites. For the structure model, Bosch chose a complete finite element modeling of the wing, the canopy is modeled by triangular shell elements and the battens by beam elements. The fluid model is the same as the one used by Breukels [START_REF] Breukels | An engineering methodology for kite design[END_REF]) and gives the forces distribution on the canopy as a function of the angle of incidence, thickness and camber of the section. The dynamic behavior of the wing is determined in an outer loop, where the tethers and the bridles are modeled. Bridle lines are modeled by linear spring-damper elements and tethers by a distance constraint. The deflection of the tethers are then neglected. For each integration step, the structure model in interaction with the fluid model gives the resultant forces applied by the wing on the bridle lines.

The results of the outer dynamic loop are set as displacement boundary conditions on the attachment point of the bridles in the structure model. This system allows the modeling of a kite in flight, for example a 8 shape-trajectory. However, it should be noticed that Bosch works is based on specific kite wings called C-shape, with only four bridle attachment points, at the tips of the wing.

All these cases are examples of partitioned Fluid-Structure Interactions where the fluid model and the structure model are segregated. An other approach to FSI problems is the monolithic coupling where fluid and structure are solved simultaneously [START_REF] Tallec | Fluid structure interaction with large structural displacements[END_REF]). Monolithic coupling is more stable and more accurate than partitioned coupling [START_REF] Michler | A monolithic approach to fluid-structure interaction[END_REF]) but also more computationally expensive. Furthermore, partitioned coupling allows the development of each model separately and is then easier to adapt to various problems. This option was therefore chosen. Nevertheless, an iterative process is necessary due to the range of deformations to deal with.

Fluid-structure coupling

This Section is devoted to the presentation of the interactions between the fluid model and the structure model. Bridles and tethers models are presented in Section 7.2.1, the data exchange between the two models is presented in Section 7.2.2.

Bridles and tether implementation

The aim of this Section is not to develop the computation of the aerodynamic forces on the tethers (see Section 4) but to describe the tether from a structural perspective. The rope used for kite applications in the beyond-the-sea® project are made of Dyneema, a lightweight and high-strength fiber. Tethers are the link between the kite and the ground, cables which allow the control of the kite. The set of bridles has also a significant role to keep the kite in shape and to transfer stresses between the kite and the tethers. As a light and very long cable, the model of the tethers triggers different challenges if the kite is modeled in static or in dynamic flights.

For example, Williams and Breukels worked on the dynamic behavior of the kite.

Williams (Williams et al. (2007c)) modeled the tethers with a series of point masses connected by inelastic links. This allows to take into account the tether drag function of the relative velocity of the wind, gravity and the aerodynamic forces on the kite, applied at the end of the tethers. However the elastic properties of the tethers are neglected.

Breukels [START_REF] Breukels | A multi-body dynamics approach to a cable simulator for kites[END_REF]) models the tether as a chain of discrete elements with a mass and an infinite stiffness, that means the elongation of the tether is neglected. The cable elements are linked with hook joints and can rotate in every direction beyond their own direction axis. Indeed, twisting has only a small effect on the flight dynamics of the kite and would force a very small time integration step. Breukels takes into account the aerodynamic drag on the tethers but also a material-based damping, which depends of the tension in the cable. The model is compared to an analytical and an experimental case of a simple pendulum and shows good agreement for both. Breukels shows that the aerodynamic drag has the most significant effects on the tether behavior.

Groot (de [START_REF] De Groot | Modeling Kite Flight Dynamics Using a Multibody Reduction Approach[END_REF]) on the contrary, chose to model the tethers as a massless straight line, neglecting the sag of the cable. The aerodynamic drag is also neglected as Groot affirmed that this approximation is valid for short tethers, up to 100m.

For the following studies, the simulation case of the kite will be static or quasi-static. Therefore, the dynamic behavior of the tethers are not modeled. In order to simplify the inputs of the structure model, the set of bridles and the tethers are modeled by straight elements, denoted by T3D2 in Abaqus™ . This type of element has no bending stiffness, no transverse shear stiffness and transfers axial force only. The compression stiffness can also be suppressed, which is the expected behavior of a line without tension. As for simplification purpose the tether is just one straight element, the aerodynamic drag is not added as a load in the structure model and is therefore neglected. These modeling choices allow to get the deformation and strain in the wing and also the tension in the tethers and the set of bridles.

For the IDK model, the complete finite element model, the tethers and the set of bridles are connected to the nodes of the leading edge structure mesh, or the trailing edge for the back tether (see Figure 7.1). For the KaB model, as the wing is modeled as an assembly of beams, the bridles are connected to rigid bodies (see Figure 7.1), which themselves are connected to the micro-beams. The length of the rigid bodies is function of the local chord length and thus the wing and its set of lines can be modeled.

Data exchange

As presented before in Section 7.1, the interaction between the fluid and structure model in our case is an exchange of deformation data and aerodynamic loading. It should be noticed that, even if the distortion of the canopy is available with the IDK model, the aerodynamic section is supposed rigid for both models [START_REF] Maison | De la problèmatique du dimensionnement de grands kites[END_REF]), that means that the set of aerodynamic coefficients do not take into account the deformations of the wing.

FSI with the Kite as a Beam structure model Simplicity was the leading principle during the development of the KaB model. For that reason, the mesh of the structure model is based on the fluid mesh (see Section 5.3).

After the computation of the fluid model, aerodynamic forces and moments are known in each point M i of the lifting line, corresponding to the tips of the micro-beams. The loads are then implemented in the structure model. It should be noticed that the initial geometry in the structure model is always the same and is only modified for each fluid computation. After convergence of the structure iteration, displacements and rotations are known for each point M i . As the kite is modeled in Abaqus™ based on its bridles and tethers, there is a rigid body motion. In order to remove this movement from the fluid iteration, the rotation and displacement of the point K, point at the quarter chord in the symmetry plane of the wing, is computed. By deducting the point K displacement from the global movement of the points M i , only the deformations of the wing are taken into account for the fluid computation.

M def i = M init i + # » U (M i ) - # » U (K) (7.1)
The rotation computed at M i is then applied to the whole horseshoe vortex. Abaqus™ gives the rotation in the form of a vector #» r , the axis of rotation and of norm η, the angle of rotation. The rotation of the horseshoe vortex is computed by following Equation 7.2:

P def = M def i + cos η # » M i P + (1 -cos η) # » M i P . #» r #» r + sin η #» r × # » M i P (7.2)
Where P is the M i corresponding point from the leading or trailing edge, or from the horseshoe vortex (see Figure 2.6 and Section 2.2).

This rigid body motion, and especially the rotation of the point K is also used after the first structure iteration for the loading of the wing. Indeed, in Abaqus™ , the loads described by force and moment vectors do not follow the deformations of the structure. For that reason, for each structure iteration k, the local aerodynamic loads are implemented on the initial geometry, non deformed, but in the converged coordinates of the structure iteration k -1. This means that a rotation corresponding to the rigid body motion of iteration k -1 is applied on the aerodynamic loads of iteration k.

FSI with the complete finite element model IDK

The interaction between the complete finite element model IDK and the nonlinear lifting line is obviously more complex since the two meshes do not even coincide. The output of the fluid model is here the distribution of effective angle of incidence α prj along the span. As a reminder, the effective angle of incidence is defined as the angle between the chord of the section and the projected effective wind, which is actually the local effective wind projected in the plane of the section. (see Section 2.2.3 and Figure 7.2).

The aerodynamic loading in the IDK model consists to apply a pressure on the canopy and forces and moments on the nodes of the leading edge. The method presented thereafter to determine the aerodynamic loading in the IDK model is an evolution of the FSI couplings presented by [START_REF] Leloup | Model approach and numerical tool development for kite performance assesment and mechanical design; application to vessels auxiliary propulsion[END_REF]) and Maison [START_REF] Maison | De la problèmatique du dimensionnement de grands kites[END_REF]). In order to estimate the aerodynamic coefficients, required for the lifting line model, and the pressure along the section, a set of 2D simulations on a Leading Edge Inflatable (LEI) kite section has been carried out with XFOIL from angle of incidence ranging from -4°to 19°. The effective angle of incidence is given by the 3D Non-Linear Lifting Line (LP3DNL, see Section 2.2) model and the pressure coefficient curve is determined from the set of 2D simulations with a cubic interpolation. As a strict LEI kite section is hard to model because of the junction between the tube and Subsequently, the pressure on the shell element of the canopy has to be evaluated. It is assumed that the pressure in the recirculation area is constant, with a value equal to the pressure at the detachment point S D (see Figure 7.4). For the first part of the canopy (X ∈ [X J ; 1] in Figure 7.3), the pressure coefficient is deduced directly from the XFOIL curve:

M A F #» x s #» z s #» t # » x prj α prj #» V prj #» V prj = #» V e -( #» V e . #» t ) #» t #» V e = # » V a + # » V ind
C p (X ∈ [X J ; 1]) = C p (Top side) -C p (Bottom side) (7.3)
For the second part of the canopy, X ∈ [X S ; X J ], as the pressure is constant in the recirculation area, the pressure coefficient is given by the following equation: Each node of the structure mesh can therefore receive a pressure coefficient value. The local pressure P on the shell element is then calculated as:

C p (X ∈ [X S ; X J ]) = C p (Top side) -C p (S D ) ( 7
P = 1 2 ρV 2 prj 1 4 (C i p + C j p + C k p + C l p ) (7.5) with # »
V prj the projected effective wind and C i p , C j p , C k p and C l p the respective pressure coefficient of each node of the quadrilateral shell element.

The last step is to compute the forces # » F LE on the leading edge node, that means the loading on the inflatable tube. XFOIL provides the pressure coefficient evolution from the sewing point S S to the detachment point S D . The remaining part of the tube is considered in the recirculation area and subjected to a constant pressure.

# » F LE = 1 2 ρV 2 prj S D S S C p (S) #» n ds + C p (S D ) S S S D #» n ds cdl (7.6)
With #» n the surface normal to the tube, s the curvilinear abscissa of point S around the tube, c and dl the local chord and the local width of the mesh respectively. As with the Kite as a Beam model, a rotation is applied on the force and moment vectors, rotation corresponding to the rigid body motion of the previos structure iteration.

After convergence of the structure iteration, the deformation of the leading edge and the trailing edge are known. The rigid body motion is here again removed from the global deformation and a new fluid mesh is computed from the distorted leading and trailing edge.

Comparison between IDK and KaB models

This Section presents a FSI computation carried out with both structure models. The wing used for the comparison tests is a kite of 5m 2 designed by beyond the sea®. This wing has a classic kite geometry, with twist and a nonlinear swept law. The case of study has been chosen from an experimental case conducted by [START_REF] Behrel | Investigation of kites for auxiliary ship propulsion: experiment setup, trials, data analysis and kite specs novel identification approach[END_REF]).

In this case, the wing was flying in a dynamic flight on a 8 shape-trajectory with a tether length of 50m. The computation case is extracted from the middle of the flight, where the wing has almost a straight trajectory. In this situation, the apparent wind seen by the wing is estimated at 20.45m/s. The traction exerted by the wing was measured and reached 737N . Apparent wind and force magnitude values give an order of magnitude for the following simulations.

As the purpose of the IDK FSI model is to predict the stress distribution in the kite for various angles of incidence in the symmetry plane, the KaB FSI model has to follow the same procedure. An angle of incidence is chosen and for each fluid iteration, the angle of incidence of the chord in the symmetry plane of the wing is set at the chosen angle. In the IDK model, based on symmetry considerations only the half of the wing is modeled. The translation along #» y and the rotation around #» x and #» z of the chord in the symmetry plane are blocked. To avoid any non-symmetrical motion, the translation along #» y and rotation around #» x and #» z of the point K in the symmetry plane of the KaB model are also blocked. Eventually, the three translations of the attachment points of the tethers on the ground are blocked (see Figure 7.5). The pressure chosen for the computation of the equivalent beam properties is 150P a, estimated from experimental results. As presented in Section 5.2, the computations are carried out in Abaqus™ , with an implicit dynamic solver applied for a quasi-static case.

Each FSI models take around 3 -5 iterations to converge as it can be seen in Figure 7.6. Each IDK structure iteration takes in average 1500s of CPU time whereas the computation time takes only of 21s for the KaB structure model, on a desktop computer (Intel E3-1220, 4 cores, 3.10GHz, 4Go RAM). Figure 7.6 shows also the necessity of a FSI kite model as the loading of the kite can increase by about 40% in some cases. IDK model. The relative variation of chord length from the initial geometry is around 5 -6% at the tips of the wing and 1 -2% in the symmetry plane. Even if the models disagree on the opening of the wing, it should be noted that the deformations which can be seen in Figure 7.8 are not very large. Finally the twist of the wing, which plays a significant part for the aerodynamic loading, is detailed in the following Section.

K UR x = UR z = 0 U y = 0 #» y #» z #» x #» z

Local results

Figures 7.9, 7.10 and 7.11 show the most significant local results on the wing. Figure 7.9 is the geometrical twist after convergence along the span of the wing, which depends of the initial geometry and of the deformations of the wing. Figure 7.9 shows a major difference between the behavior of the two models. The local twist of the KaB model is almost constant at -3°in the center of the wing (value of the initial twist) and increases up to 11.5°at the tips. After convergence, the twist of the KaB wing varies only at the tips. On the contrary, the local twist of the IDK model decreases almost linearly from -3°in the center to -32°at the tips. Figure 7.10 shows the projected effective angle of incidence and the projected wind velocity on the local sections after convergence. As defined in the LP3DNL model, a negative twist induces a positive angle of incidence. The effective angle of incidence is directly linked to the twist of the wing and depends also on the local induced wind (see Section 7.2.2 or Section 2.2 for the definition of the effective angle of incidence). The differences shown in Figure 7.9 are therefore propagated in Figure 7.10. The effective angle of the KaB model is almost constant around 1°over most of the span of the wing and decreases until -8°at the tips. The effective angle of the IDK model is in the range [6°; 9°] in the middle of the wing and increases sharply up to 22°and eventually 35°at the tips of the wing.

The projected wind velocities of the two models are in the same range of values, the skip in the IDK effective angle can be found again in the effective wind curve around the abscissa -1m and 1m, even if it is less pronounced.

Eventually, Figure 7.11 shows the magnitude of the local loading per unit length along the span of the wing. The effect of the differences in the local twist can be found here also. The local aerodynamic loading per unit length is overall larger for the IDK model, as expected with the global aerodynamic loading presented in Section 7.3.1 and also with the effective angle of incidence of

Discussion

The differences between the two models observed in Chapter 6 are still present in Fluid-Structure Interaction results, the set of bridles is not enough to smooth out the gap between the two structure models. The two major geometrical differences are the wingsweep and the twist at the tips of the wing which can be seen in Figures 7.7 and 7.9 respectively. The most significant difference is the resulting twist of the wing which increases for the KaB model and decreases to large angles for the IDK model. Such large angles are quite unexpected however, in reality, the kite designed by beyond-the-sea® has a pulley between the front and back tether which is not modeled here. This pulley, fastened near the tips of the wing, transfers part of the stress between the front and back tethers and would come into play in our case since the equivalent bridle in the simulation is not loaded. Secondly for the pressure loading, the IDK model works with a set of pressure coefficient distribution computed with XFOIL for a limited range of incidence. If the angle of incidence of the section is outside of this range, the lower or upper limit is taken, which limits the influence of the angle of incidence when the section is stalled and may prevent to come back to lower angles of incidence. The upper limit of the computed range of incidence for the previous simulations was 19°, beginning of the stalling range of the section. It will be certainly useful to have a wider range of angles of incidence for further simulations.

A first source of differences between the two models may be the way the bridles are attached to the wing. For the IDK model, the tethers are fastened to a beam element which models the inflatable leading edge. The volume of the tubes are not represented.

For the KaB model, the tethers are linked to the wing beams by a rigid body and of the length of the complete LEI kite section. Figure 7.12 shows that the tethers are not attached exactly at the same position. However, as the fluid computation was done with a constant angle of incidence, this does not seem to be the major source of difference between the two models. Figure 7.13 shows the local twist for the initial geometry, after the first fluid-structure iteration and after convergence for the IDK and the KaB models. After only one iteration the twist is already very close to the converged twist. Even if the lifting line model has not been validated for high twist angles, after only one iteration with an equivalent aerodynamic loading, the models already differ significantly. The same phenomena can be found for the other deformations (wingsweep or opening of the wing). Section 6.2 has shown the gap between the two models for the modeling of the twist, as the torsion stiffness of the KaB model is significantly higher. Figure 7.13 and the previous results of Section 7.3 show the consequence of this overestimation.

The second major difference is the wingsweep of the wing, the displacement along the #» x axis. The two models deform in the same direction, towards the reduction of the initial wingsweep. However the distortion of the KaB model is definitely bigger. This is still consistent with the results of Section 6. model is 150P a, the stiffness of the IDK model along the #» x axis is then also bigger than the KaB one.

Section 6.2 also showed a large difference in stiffness along the #» y and #» z axis but these stiffnesses do not seem to play a critical role in our case.

Mesh influence on the KaB model

For the previous simulations, a mesh with 30 sections per half-span (60 sections on the whole wing) was chosen as this discretization is a good compromise between robustness and accuracy for the 3D Non-Linear Lifting Line model (see Section 2.3.2). However, as the density of the mesh has an influence on the aerodynamic loading (see Figure 7.14), the influence of the mesh on the Kite as Beam fluid-structure interaction has to be investigated. All following results are for the previous kite wing at 7°of incidence.

Figure 7.14 shows in red the aerodynamic loading on the wing as a function of the number of sections for the converged result and the first iteration (with an undeformed wing) in black. The relative difference is computed by considering the finest mesh of 150 sections as reference. As the structural problem stays relatively simple even when refining the mesh, the computation time does not increase significantly.

As it can be seen in Figure 7.14, the estimation difference on the aerodynamic loading for the first iteration is quite small, less than 10% between 20 and 150 sections over the span of the wing. However, this has a large influence on the converged result.

There is almost a 50% gap between the 20 sections and the 150 sections mesh for the aerodynamic loading. It should be noticed that increasing the number of sections increases the drag of the wing (18% difference between 20 and 150 sections) and cut 105 in half the lift of the wing. To simplify the coupling between the two models, it was chosen to work with the same mesh density for the fluid model and the structure model. But Figure 7.14 shows that the influence of discretization should be checked for the KaB structure model.

Figure 7.15 shows the local twist of the sections and the local aerodynamic per unit length over the span of the wing for 20, 60 and 150 sections. The number of sections obviously changes the converged results, the 20 sections wing is overall more loaded, but the trends over the span are similar. As the IDK model was discretized with 20 sections over the span and as the number of sections has an influence at least on the aerodynamic loading, it was chosen to work with this mesh density for the following studies.

7.4.2 Influence of the work pressure on the KaB model

Previous simulations have all been carried out with the same work pressure, 150P a, estimated from experimental results. However in the IDK model, the pressure is recalculated after each fluid iteration. To investigate the influence of the work pressure on the results, the pressure over the span has been extracted from the converged results of the IDK FSI. Each elementary cell has then a distinct pressure between [250; 300P a], averaged from the pressures of the equivalent wing section in the IDK model. The global results obtained are compared in Table 7.2 with the results from the IDK FSI and the wing with 20 sections and a work pressure of 150P a, to serve as a reference.

The converged aerodynamic results are very close and the only significant difference is the deformation of the wing. The local results and the deformed geometries can be found in Figure 7.17. As expected, the KaB model with increased pressures has a global increased stiffness. The wingsweep of the wing at the tips is decreased as well as the opening of the wing. The work pressure has then a small influence, not on the aerodynamic loading (the local loads are quasi identical in Figure 7.17 Aerodynamic loading (N)

Possibilities of model calibration

As the KaB model is really away from the IDK model, a last computation case has been done to evaluate its calibration possibilities to improve the results. Since the twist and the wingsweep are the two largest differences between the two models, it was chosen to decrease the torsional stiffness of the equivalent beams GJ by a coefficient 0.5 and to increase the in-plane bending stiffness EI 3 by a coefficient 2.0. It should be noticed that the coefficients were constant for each equivalent beam while geometrical approximations between the two models are the largest at the tips of the wing. The results presented in Table 7.3 can certainly be improved by adjusting more precisely the equivalent beams stiffnesses. The local results can be seen in Figure 7.18.

The global aerodynamic loading increases by around 50% and the wingsweep decreases from about 15%. However the global opening of the wing also increases significantly, which is coherent with the modification of the stiffnesses. Indeed, decreasing the torsional stiffness, increases significantly the twist at the tips of the wing, which increases the local aerodynamic loading and thus the opening of the wing.

The results can therefore be improved by the use of coefficients to calibrate the KaB model. However this implies the computation of at least the comparison cases presented in Section 6.1 for each new geometries to adjust precisely the coefficients. These simulations cases are not very time consuming, considering the structure model is already build, but this reduces the ease of use of the Kite as a Beam model.

IDK model

Conclusion

The aim of this study was to develop a fast and reliable modeling approach, able to predict the structural deformations, the tension in the tethers and the aerodynamic performances of a leading edge inflatable kite according to the wind direction and the wing kinematic. Two models have been developed and also, the interface needed to a fluid-structure interaction simulation loop. Each model has been compared to a more complex and sophisticated one, for validation purposes.

A 3D Non-Linear Lifting Line model has been developed. This model is based on the lifting line of Prandtl, extended to deal with non straight kite wings. The model uses non linear aerodynamic coefficients which can easily be extracted from experiment or simulation data. The finite wing and its wake are represented by a set of horseshoe vortices which consists of six vortex segments. For each horseshoe vortex, the bound vortex is orthogonal to the two adjacent trailing vortices, which extend parallel to the chord over one chord length. The second part of the trailing vortices extends parallel to the global relative free stream direction. An algorithm was implemented to calculate the circulation of each horseshoe vortex and then the lift, drag and moment of the wing. A damping factor is necessary for the good convergence of the algorithm. Numerical studies showed that a wake length equivalent to 20 times the chord in the symmetry plane of the wing and a discretization between 30 and 50 sections per halfwing provides converged results in the numerical accuracy range of the model. Finally, the model showed perfect agreement with the lifting line theory for classic wing geometries.

To further validate the 3D Non-Linear Lifting Line, RANSE simulations were performed. At first, 2D simulations were carried out to compute the aerodynamic coefficients used in the lifting line model and also to estimate the numerical accuracy of the RANSE simulations. 3D simulations have been carried out on three geometries to test the limits of the lifting line model. For the geometry close to a leading edge inflatable kite, lift coefficient estimations have less than 5% difference over a large range of angle of incidence. The estimation of the drag coefficient is also satisfying even if the gap is a bit more significant for angles of incidence greater than 10°. Simulations with the wing with a sideslip angle showed that the lifting line model is still able to predict correctly to evolution of the aerodynamic coefficients. Local results were also computed on the 3D RANSE simulations. The lifting line model is able to predict the position of minimum and maximum loading, for a wing with various angles of incidence and a sideslip angle in the range of [0°; 15°].

The model was also introduced in a new equilibrium iterative procedure for the kite, able to calculate the effective lift-to-drag ratio of a kite along a circular flight path, taking into account its turn rate. A parametric study was conducted to analyze the evolution of the aerodynamic coefficients and the flight characteristics according to the tether length or the radius of the trajectory. The results showed the influence of the turn rate on the aerodynamic loading of the kite and the capacity of the lifting line model to take this parameter into account.

Simultaneously, a structure model was developed. The core idea of the Kite as a Beam model is to approximate a Leading Edge Inflatable kite by an assembly of equivalent beam elements. Over the span, the kite is divided into several elementary cells, each one composed of a portion of the inflatable leading edge, two inflatable battens and the corresponding canopy. A cell is made up of beam elements for the battens and the leading edge and of shell elements for the canopy. If the wing has a leech line, it is modeled with truss elements at the trailing edge. First of all, the cell is put under pressure with its four corner clamped. The cell is then subjected to different linear displacement perturbations applied on one side while the other side is clamped. The equivalent beam stiffnesses are calculated from the reaction forces and moments computed after each linear perturbation. Eventually, the equivalent beam are re-discretized in order to match with the lifting line mesh.

The Kite as a Beam model was then compared to a complete kite Finite Element model. The structure models were put under a constant pressure with equivalent boundary conditions. The same elementary displacements were then applied, a translation along one of the three global axes or a rotation of the wing tip. The reaction forces and moments were compared for two different wings, one modeling a kite designed for the beyond the sea® project and an academic wing taken cylindrical with a constant chord. Results showed that the Kite as a Beam model largely overestimates the torsion for kite geometries, where large geometrical approximations are made when building the elementary cell. The Kite as a Beam model was also not able to follow the evolution trend of the complete Finite Element model for the opening of the kite, even if it was in the order of magnitude of reaction forces.

Nevertheless, the Kite as a Beam model was coupled with the 3D Non-Linear Lifting Line model to set the basis of a FSI design loop. The fluid model computes the local aerodynamic load on the wing, which are applied on the structure model. After convergence of the Kite as a Beam model, the deformations are computed by subtracting the rigid body motion of the wing and are taken as input for the lifting line model. The same principle was used for the complete Finite Element model, except that the wing is loaded through pressure on the canopy. The gain in computation time is really significant as the Kite as a Beam model takes only a few minutes to compute while the complete Finite Element model needs one or two hours. The results showed the influence of the mesh density on the converged results, especially for the aerodynamic loading, which should be kept in mind if the model is used to provide quantitative values. The comparison between the two fluid-structure models pointed out their behavior differences, but also the fact that these differences could be reduced by calibrating the torsional stiffness and the in-plane bending stiffness of the Kite as a Beam model.

Perspectives

The fluid-structure simulations showed the influence of the mesh density on the Kite as a Beam results. To simplify the coupling between the two models, it was chosen to work with the same mesh density for the fluid model and the structure model, but a mesh convergence study should also be done for the structure model only, to estimate the impact of the number of beam elements without the influence of the fluid model. The model could also take into account the pressure resulting from the fluid iteration, with a linear interpolation between two or three equivalent beam properties computed at different pressure. However, the purpose of the Kite as a Beam model is to be able to test rapidly different wing geometries, shapes and sizes, for various wind conditions. The results are more interesting from a relative point of view than as for quantitative values. From this perspective, the calibrated Kite as a Beam model should be compared to the complete Finite Element model for another geometry (with more wingsweep for example) or another angle of incidence, to check if the calibrated model is able to follow the evolution trend of the Finite Element model. If the results of the Kite as a Beam FSI model are not enough for a first estimation in a design phase, then another model should be developed with less strong assumptions, for example for the tangent stiffness of the elementary cell.

For the fluid part, the 3D Non-Linear Lifting Line model has been checked for various angles of incidence and of sideslip, but only in the steady range. Other simulations could be carried out in the unsteady range or with high twist angles, to further validate the results of the lifting line. Another simulation could also be done with a LEI kite profile, to check its influence on the lifting line estimation. In the same way, if the fluid model is able to take into account the turn rate of the kite, no validation of the global and local aerodynamic loading of the wing with turn rate has been provided. Lastly, the lifting line model uses only one set of aerodynamic coefficients. One could add the option in the program to use another set of aerodynamic coefficients, to take into account an evolution of the aerodynamic section over the span of the wing, or the section of the battens.

Finally for the global fluid-structure interaction, if the Kite as Beam model is very different from the complete Finite Element model, one can not assess on the validity of either models. Indeed, the twist of the complete Finite Element model was really significant and the deformations observed during real experiments were not as severe.

To further validate the Fluid-Structure Interactions models, accurate information on the deformations of the wing would be greatly useful. For example, Augier [START_REF] Augier | Etudes expérimentales de l'interaction fluide-structure sur surface souple: application aux voiles de bateaux[END_REF]) used a camera to obtain the deformations of a boat sail and validates a model of fluid-structure coupling or Breukels [START_REF] Breukels | An engineering methodology for kite design[END_REF]) carried out experiments with a LEI wing in wind tunnel where photogrammetry was used to get the shape of the canopy. The second research axis is the modeling of the FSI directly for a kite in dynamic flight. At the moment, the kinematic of the kite, its velocity and its turn rate, are inputs of the fluid model. Modeling the kite in dynamic flight would allow to estimate the traction of the kite according to its geometry, its trajectory and the apparent wind all along its trajectory. One could then better compare these results with the real experiments carried out by [START_REF] Behrel | Investigation of kites for auxiliary ship propulsion: experiment setup, trials, data analysis and kite specs novel identification approach[END_REF]).

Section coefficients

To finish the description of the wing, one needs to define the aerodynamic section, which is described by three curves in the lifting line model, the lift coefficient, the drag coefficient and the moment coefficient curves. It was chosen to work with parametric curves, to take into account the non-linearity of the lift coefficient. The parameters are chosen to facilitate the reproduction of curves obtained with experiments or CFD simulations, by just reading the graph.

Lift coefficient

Part I, II and V of the lift polar curve are linear (see Figure A.1), while part III and IV are described by a polynomial of degree 3 and 2 respectively. The different curve parts are tangent in (α 2 , C l 2 ) and in (α 3 , C l 3 ). The points (α i , C l i ) divide the different parts of the curve which will pass through each point. (α 3 , C l 3 ) is the maximal value of the lift polar curve. To describe the lift polar curve, one only needs these four points and two coefficients for the slope of the linear parts I and V.

1. 1

 1 Left, painted ceramic disc representing a boat with two mast Carter (2006). Right, De Vlijt windmill in the Netherlands . . . . . . . . . . . 1.2 Left, kite from the Bali Kite Festival. Right, kite from the beyond the sea® project . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Left, catamaran Race For Water, project which fight against plastic pollution in the oceans. Right, catamaran Energy Observer, project seeking to develop solutions for the energy transition. . . . . . . . . . . 2.1 Induced velocity by a vortex filament of strength Γ Anderson (2011) . 2.2 Prandtl's lifting line model (Anderson (2011)) with L the lift and D i the induced drag on the wing . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Definition of the lifting line . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Lift polar curve. Parts I, II and V are linear while parts III and IV are polynomial (3rd and 2nd degree respectively) . . . . . . . . . . . . . . 2.5 Example of a low discretized lifting line model. Red vectors are proportional to the local aerodynamic resulting forces. To improve the readability local torques were not represented. . . . . . . . . . . . . . . 2.6 Section geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.7 Local effective wind projected in the section plane . . . . . . . . . . . . 2.8 Initial and converged circulations along the span of the wing . . . . . . 2.9 |C X (N )-C X (10 3 )| |C X (10 3) | Relative error on the aerodynamic coefficients with respect to the wake length nondimensionalized by the root chord length. 2.10 (a) Averaged computation time with respect to the number of sections per half-wing. (b) Number of iterations vs. the number of sections per half-wing. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.11 Evolution of the relative variation of the aerodynamic coefficients with respect to the number of sections per half-wing. The relative variation is estimated as |C X (N )-C X (N -1)| 10|C X (N )|+1 , using the results of two consecutive cases, N -1 and N . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.12 Wing with an elliptic plan form and an elliptic circulation distribution . 2.13 (a) Circulation along the span of the wing. (b) Lift coefficient with respect to the angle of incidence . . . . . . . . . . . . . . . . . . . . . . 2.14 (a) Spanwise circulation for a wing at 5°of incidence (b) Lift coefficient with respect to the angle of incidence. Wing with a maximal twist of 5°in red, or -5°in blue. . . . . . . . . . . . . . . . . . . . . . . xi 3.1 2D mesh of the computational domain . . . . . . . . . . . . . . . . . . 3.2 (a) View of the mesh around the profile section. View of the near wall mesh in the vicinity of (b) the leading edge and (c) the trailing edge. . . 3.3 The 3D computational domain . . . . . . . . . . . . . . . . . . . . . . . 3.4 (a) Slices in the wing geometry, (b) 3D sections and mesh refinements in the wake, around the tip vortex and the slices . . . . . . . . . . . . . 3.5 Total force per unit length on each slice function of the slice thickness, (a) at the wing root, (b) in the middle of the half-wing, (c) at the tip. In gray, the estimated converged value. The blue triangle is the chosen thickness. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 2D polar curves obtained via STAR-CCM+® (filled squares), compared with experimental data Abbott & Doenhoff (1959) (filled triangles), and fitted with parametric polar curves used for the lifting line method (plain line). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.7 Convergence history of the lift coefficient for variations (a) of the computational domain size, (b) of the number of cells in the mesh, (c) of the number of layers of the near wall mesh and of the turbulence model. 3.8 Convergence history of the drag coefficient for variations (a) of the computational domain size, (b) of the number of cells in the mesh, (c) of the number of layers of the near wall mesh and of the turbulence model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.9 Un-twisted and un-swept geometry (case 1), (a) side view, (b) front view, (c) developed view . . . . . . . . . . . . . . . . . . . . . . . . . . 3.10 3D aerodynamic coefficients as a function of angle of incidence, obtained via STAR-CCM+® (filled squares) and the lifting line method (plain line) on the wing (case 1). Error bars are the estimated numerical accuracy. No angle of sideslip. . . . . . . . . . . . . . . . . . . . . . . . 3.11 3D aerodynamic coefficient as a function of sideslip angle obtained via STAR-CCM+® (filled squares) and the lifting line method (plain line) on the wing (case 1) at 2°of incidence. Error bars are the estimated numerical accuracy. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.12 Magnitude of the local aerodynamic force per unit length on the case 1 wing obtained via STAR-CCM+® (filled squares) and the lifting line method (filled triangles). (a) Wing at 2°of incidence and 0°of sideslip. (b) Wing at 12°of incidence, 0°of sideslip. (c) Wing at 2°of incidence and 7.5°of sideslip. (d) Wing at 2°of incidence and 15°of sideslip. . . 3.13 Un-twisted and linearly swept geometry (case 2), (a) side view, (b) front view, (c) developed view . . . . . . . . . . . . . . . . . . . . . . . 3.14 3D aerodynamic coefficients as a function of angle of incidence, obtained via STAR-CCM+® (filled squares) and the lifting line method (plain line) on the wing (case 2). Error bars are the estimated numerical accuracy. No angle of sideslip. . . . . . . . . . . . . . . . . . . . . . . . xii 3.15 Magnitude of the local aerodynamic force per unit length on the case 2 wing obtained via STAR-CCM+® (filled squares) and the lifting line method (filled triangles). No angle of sideslip. (a) Wing at -4°of incidence. (b) Wing at 0°of incidence. (c) Wing at 10°of incidence. (d) Wing at 15°of incidence. . . . . . . . . . . . . . . . . . . . . . . . . 41 3.16 Magnitude of the local aerodynamic force per unit length on the swept wing (case 2) obtained via STAR-CCM+® (filled squares) and the lifting line method (filled triangles). (a) Wing at 5°of incidence and 0°o f sideslip. (b) Wing at 5°of incidence and 7.5°of sideslip. . . . . . . 42 3.17 Linearly twisted and non-linearly swept geometry (case 3), (a) side view, (b) front view, (c) developed view . . . . . . . . . . . . . . . . . . 43 3.18 3D aerodynamic coefficients as a function of angle of incidence or angle of sideslip (at 5°of incidence), obtained via STAR-CCM+® (filled squares) and the lifting line method (plain line) on the kite wing (case 3). Error bars are the estimated numerical accuracy. . . . . . . . . . . 44 3.19 Magnitude of the local aerodynamic force per unit length on the kite wing (case 3) obtained via STAR-CCM+® (filled squares) and the lifting line method (filled triangles). No angle of sideslip. (a) Wing at -4°o f incidence. (b) Wing at 0°of incidence. (c) Wing at 10°of incidence. (d) Wing at 15°of incidence. . . . . . . . . . . . . . . . . . . . . . . . . 45 3.20 Magnitude of the local aerodynamic force per unit length on the kite wing (case 3) obtained via STAR-CCM+® (filled squares) and the lifting line method (filled triangles). (a) Wing at 5°of incidence and 0°o f sideslip. (b) Wing at 5°of incidence and 7.5°of sideslip. (c) Wing at 5°of incidence and 15°of sideslip. . . . . . . . . . . . . . . . . . . . 46 4.1 Circular path of the kite, with V RW the wind, x V K the kite velocity direction, V a the apparent wind seen by the kite, O the anchor point and R the radius of the trajectory. . . . . . . . . . . . . . . . . . . . . . 50 4.2 (a) Non-dimensional norm of the local aerodynamic resultant along the kite as obtained from the equilibrium procedure, with and without taking into account the turn rate influence. The local norms are divided by the average value of the local norms in the case of the kite in pure translation motion. (b) General flight specs. . . . . . . . . . . . . . . . 52 4.3 Evolution of (a) the lift-to-drag ratio, (b) the lift coefficient and (c) the drag coefficient depending on the radius of the trajectory divided by the tether length, and for different lengths of tether from 25m to 200m. 53 4.4 Evolution of (a) the roll angle, of (b) the kite velocity and of (c) its turn rate depending on the radius of the trajectory divided by the tether length, for different lengths of tether from 25m to 200m. . . . . . . . . 54 5.1 Denomination of the main kite elements . . . . . . . . . . . . . . . . . 57 5.2 Buckling of a LEI kite during flight. Picture taken during experimental tests of Behrel Behrel (2017) . . . . . . . . . . . . . . . . . . . . . . . . 59 xiii 5.3 Real profile of a LEI kite. The black cross denotes the position of the sewing between the inflatable tube and the canopy. The blue curve represents the canopy in the IDK model. . . . . . . . . . . . . . . . . . 5.4 Example of a LEI kite according to the IDK model . . . . . . . . . . . 5.5 Kite as a Beam concept and division of the kite in several elementary cells . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.6 Elementary cell geometry . . . . . . . . . . . . . . . . . . . . . . . . . 5.7 Geometry of the elementary cell . . . . . . . . . . . . . . . . . . . . . . 5.8 Shape of the canopy under homogeneous pressure loading with contour value representing the out-of-plane displacement along - → z c (in meter) 5.9 Case (β): out-of-plane shear along z c . The contour value represents the displacement along z c . . . . . . . . . . . . . . . . . . . . . . . . . . 5.10 Micro-beams and equivalent beams with the kite respective geometry . 5.11 Determination of the centroid of a micro-beam . . . . . . . . . . . . . . 5.12 Relative error on the (α) traction case with respect to the number of beam elements for the largest cell of the kite. . . . . . . . . . . . . . . . 5.13 Convergence study of the D/CS ratio . . . . . . . . . . . . . . . . . . . 6.1 Boundary conditions in displacements (U) and angular velocity (VR) of the kite under pressure load . . . . . . . . . . . . . . . . . . . . . . . 6.2 Details of the notations used for the distorted geometry equivalence between KaB and IDK . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3 Description of the different comparison cases . . . . . . . . . . . . . . 6.4 Complex Finite Element model of the LEI kite. . . . . . . . . . . . . . 6.5 Reaction forces and moment of the kite wing for the different comparison cases. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.6 Complex Finite Element model of the cylindrical wing . . . . . . . . . 6.7 Reaction forces and moment of the cylindrical wing for the different comparison cases. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.8 Deformation of one batten during the first pressure step (in red) and after the displacement step (in green). Case at 200P a and displacement along the X axis. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.9 Initial profile of the canopy . . . . . . . . . . . . . . . . . . . . . . . . . 6.10 Shape of the canopy for a coefficient R LA in range [0.25; 0.5], C MX = 0.065 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.11 Shape of the canopy for a coefficient C M X in range [0.001; 0.1], R LA = 0.3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.12 Evolution of the equivalent beam properties for position of maximum camber R LA ∈ [0.25; 0.5] for a range of pressure [100; 500P a], C MX = 0.065 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.13 Evolution of the equivalent beam properties for different camber C MX ∈ [0.001; 0.1] for a range of pressure [100; 500P a], R LA = 0.3 . . . . . . . 6.14 Equivalent beam frame and micro-beam frame from the Kite as a Beam model of a kite geometry. . . . . . . . . . . . . . . . . . . . . . . . . . . xiv 7.1 Set of bridles, tethers and connection to the wing for the IDK model (left) and the KaB model (right

  (a) Magnitude of the aerodynamic loading as a function of the number of sections on the wing in the range[20; 150] for the FSI converged result in red and for the first iteration (with an undeformed wing) in black. (b) Relative difference on the aerodynamic loading with respect to the number of sections on the full span of the wing for the converged result in red and the first iteration in black. . . . 106 7.15 Wing at 7°on incidence after convergence with 20, 60 or 150 sections over the span,(a) local twist of the aerodynamic sections, (b) magnitude of the local aerodynamic loading per unit length . . . . . . . . . . 107 xv 7.16 Geometries of the wing at 7°of incidence after convergence with 20, 60 or 150 sections over the span. (a) Solid line represents the lifting line coordinates. (b) Solid line represents the leading edge, dashed line the trailing edge. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.17 Wing at 7°of incidence after convergence with the IDK model (in red), the KaB model with a work pressure of 150P a (in blue), and with work pressures extracted from the IDK model (in green). (a) Local twist of the aerodynamic sections. (b) Magnitude of the local aerodynamic loading per unit length. (c) Solid line represent the leading edge, dashed line the trailing edge. (d) Deformed lifting line. . . . . . . 7.18 Wing at 7°of incidence after convergence with the IDK model (in red), the KaB model with a work pressure of 150P a (in blue), with calibrated properties (in green). (a) Local twist of the aerodynamic sections. (b) Magnitude of the local aerodynamic loading per unit length. (c) Solid line represent the leading edge, dashed line the trailing edge. (d) Deformed lifting line. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.1 Description of the lift polar curve . . . . . . . . . . . . . . . . . . . . . A.2 Description of the drag polar curve . . . . . . . . . . . . . . . . . . . . A.3 Description of the moment polar curve . . . . . . . . . . . . . . . . . . xvi List of Tables 2.1 Number of iterations for each of the initial circulations. . . . . . . . . . 2.2 Estimation of the confidence interval for a number of sections n between [30; 100], α is the angle of incidence of the wing . . . . . . . . . 3.1 Estimation of the confidence intervals, coming from the variations of the domain size, of the mesh and of the turbulence model. C stands for coefficient. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Overview of the different validation cases . . . . . . . . . . . . . . . . 5.1 Boundary conditions on the beam B R in displacements (U ) and rotations (U R) for the load cases (α)-(φ), components expressed in the frame ( - → x c , -→ y c , -→ z c ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Computation time (total CPU time) of all the linear cases with respect to the beam numbers for the largest cell of the kite. . . . . . . . . . . . 6.1 Boundary conditions in displacements (U), velocity (V ) and angular velocity (VR) on the IDK model for the different comparison cases . . 6.2 Boundary conditions in displacements (U), rotation (UR) and angular velocity (VR) on the KaB model for the different comparison cases . . 6.3 Average CPU time for a comparison case, the inflation step is included in the IDK simulation (usual Desktop computer) . . . . . . . . . . . . 6.4 Gap (in %) between the IDK and the KaB models for the kite geometry displacement. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1

 1 Figure 1.1 -Left, painted ceramic disc representing a boat with two mast Carter (2006). Right, De Vlijt windmill in the Netherlands

Figure 1

 1 Figure 1.2 -Left, kite from the Bali Kite Festival. Right, kite from the beyond the sea® project

Figure 1

 1 Figure 1.3 -Left, catamaran Race For Water, project which fight against plastic pollution in the oceans. Right, catamaran Energy Observer, project seeking to develop solutions for the energy transition.

Figure 2

 2 Figure 2.1 -Induced velocity by a vortex filament of strength Γ Anderson (2011)

Figure 2

 2 Figure 2.2 -Prandtl's lifting line model (Anderson (2011)) with L the lift and D i the induced drag on the wing

  Figure 2.4 shows an example of the division of the lift polar curve and the different parameters needed for the lifting line algorithm. The details on the definition of polar curves are given in Appendix 7.4.3.

Figure 2 . 4 -

 24 Figure 2.4 -Lift polar curve. Parts I, II and V are linear while parts III and IV are polynomial (3rd and 2nd degree respectively)

Figure 2

 2 Figure 2.5 -Example of a low discretized lifting line model. Red vectors are proportional to the local aerodynamic resulting forces. To improve the readability local torques were not represented.

  (w max , w * (1 + e w )) else w = max(w min , w * e w )endif(2.11)

  theorem (Eq.1.23) Weighted elliptic distribution (Eq.1.24) Mean value 2D (Eq.1.25) Converged Circulation Curves

Figure 2 . 8 -

 28 Figure 2.8 -Initial and converged circulations along the span of the wing

Figure 2

 2 Figure 2.10 -(a) Averaged computation time with respect to the number of sections per half-wing. (b) Number of iterations vs. the number of sections per half-wing.

Figure 2

 2 Figure 2.13 -(a) Circulation along the span of the wing. (b) Lift coefficient with respect to the angle of incidence

Figure 2 . 14 -

 214 Figure 2.14 -(a) Spanwise circulation for a wing at 5°of incidence (b) Lift coefficient with respect to the angle of incidence. Wing with a maximal twist of 5°in red, or -5°in blue.

  Figure 3.1 -2D mesh of the computational domain

  Figure 3.2 -(a) View of the mesh around the profile section. View of the near wall mesh in the vicinity of (b) the leading edge and (c) the trailing edge.

  Figure 3.3 -The 3D computational domain

  Figure 3.4 -(a) Slices in the wing geometry, (b) 3D sections and mesh refinements in the wake, around the tip vortex and the slices

Figure 3

 3 Figure 3.5 -Total force per unit length on each slice function of the slice thickness, (a) at the wing root, (b) in the middle of the half-wing, (c) at the tip. In gray, the estimated converged value. The blue triangle is the chosen thickness.

  Figure 3.6 -2D polar curves obtained via STAR-CCM+® (filled squares), compared with experimental data Abbott & Doenhoff (1959) (filled triangles), and fitted with parametric polar curves used for the lifting line method (plain line).

Figure 3 Figure 3

 33 Figure 3.7 -Convergence history of the lift coefficient for variations (a) of the computational domain size, (b) of the number of cells in the mesh, (c) of the number of layers of the near wall mesh and of the turbulence model.

Figure 3 Figure 3

 33 Figure 3.9 -Un-twisted and un-swept geometry (case 1), (a) side view, (b) front view, (c) developed view

Figure 3

 3 Figure 3.11 -3D aerodynamic coefficient as a function of sideslip angle obtained via STAR-CCM+® (filled squares) and the lifting line method (plain line) on the wing (case 1) at 2°o f incidence. Error bars are the estimated numerical accuracy.

  α = 2°β = 15°F igure 3.12 -Magnitude of the local aerodynamic force per unit length on the case 1 wing obtained via STAR-CCM+® (filled squares) and the lifting line method (filled triangles). (a) Wing at 2°of incidence and 0°of sideslip. (b) Wing at 12°of incidence, 0°of sideslip. (c) Wing at 2°of incidence and 7.5°of sideslip. (d) Wing at 2°of incidence and 15°of sideslip.

Figure 3

 3 Figure 3.13 -Un-twisted and linearly swept geometry (case 2), (a) side view, (b) front view, (c) developed view

Figures 3 .Figure 3 . 14 -

 3314 Figure 3.14 -3D aerodynamic coefficients as a function of angle of incidence, obtained via STAR-CCM+® (filled squares) and the lifting line method (plain line) on the wing (case 2). Error bars are the estimated numerical accuracy. No angle of sideslip.

  α = 15°β = 0°F igure 3.15 -Magnitude of the local aerodynamic force per unit length on the case 2 wing obtained via STAR-CCM+® (filled squares) and the lifting line method (filled triangles). No angle of sideslip. (a) Wing at -4°of incidence. (b) Wing at 0°of incidence. (c) Wing at 10°of incidence. (d) Wing at 15°of incidence.

  α = 5°β = 7.5°F igure 3.16 -Magnitude of the local aerodynamic force per unit length on the swept wing (case 2) obtained via STAR-CCM+® (filled squares) and the lifting line method (filled triangles). (a) Wing at 5°of incidence and 0°of sideslip. (b) Wing at 5°of incidence and 7.5°of sideslip.

Figure 3

 3 Figure 3.17 -Linearly twisted and non-linearly swept geometry (case 3), (a) side view, (b) front view, (c) developed view

Figure 3 .

 3 Figure3.18 -3D aerodynamic coefficients as a function of angle of incidence or angle of sideslip (at 5°of incidence), obtained via STAR-CCM+® (filled squares) and the lifting line method (plain line) on the kite wing (case 3). Error bars are the estimated numerical accuracy.

  [START_REF] Dadd | Determination of kite forces using three-dimensionnal flight trajectories for ship propulsion[END_REF]) used the zero-mass model to predict the kite forces evolution along an eight shaped trajectory. In this model, the weights of the kite and of the tethers are neglected. This assumption is discussed in[START_REF] Dadd | Comparison of Two Kite Force Model with Experiment[END_REF]). The results of the zero mass model are compared with a concentrated mass model, where the kite is modeled as a point mass at the end of the tether. Dadd showed that the mass of the kite and the tether can be neglected if the weight is small enough in comparison with aerodynamic forces. But he assumed that the aerodynamic coefficients are constant, as well as the angle of attack of the kite. The tether is considered straight, inextensible and dragless. Eventually, Dadd shows that the zero mass model is able to give a correct first approximation of the kite velocity in dynamic flight.The hypothesis of a straight line for the modeling of the kite is widespread (see[START_REF] Alexander | Kite equilibrium and bridle length[END_REF];[START_REF] Terink | Flight Dynamics and Stability of a Tethered Inflatable Kiteplane[END_REF]).[START_REF] Stevenson | Traction kite testing and aerodynamics[END_REF]) studied the stability of the kite during flight in a plane in order to find equilibrium points. Stevenson compared experimental measures and kite models with straight lines or curved ones. The modeling of the tethers has only a small impact on the results and Stevenson recommended the use of straight lines, even when there are several lines. Furthermore, according to[START_REF] De Groot | Modeling Kite Flight Dynamics Using a Multibody Reduction Approach[END_REF]), the tethers can be assumed straight for lengths inferior to 150m.

Figure 4 . 1 -

 41 Figure 4.1 -Circular path of the kite, with V RW the wind, x V K the kite velocity direction, V a the apparent wind seen by the kite, O the anchor point and R the radius of the trajectory.

Figure 4

 4 Figure 4.2 -(a) Non-dimensional norm of the local aerodynamic resultant along the kite as obtained from the equilibrium procedure, with and without taking into account the turn rate influence. The local norms are divided by the average value of the local norms in the case of the kite in pure translation motion. (b) General flight specs.

Figure 4 . 3 -

 43 Figure 4.3 -Evolution of (a) the lift-to-drag ratio, (b) the lift coefficient and (c) the drag coefficient depending on the radius of the trajectory divided by the tether length, and for different lengths of tether from 25m to 200m.

Figure

  Figure4.3 (a) shows the evolution of the lift-to-drag ratio depending on the relative radius, for different lengths of the tether between 25m and 200m. The lift-to-drag ratio is quasi-constant for longer tether lengths, in contrast it varies substantially for the smaller tether lengths. The same kind of variations can be found for the coefficients of lift and of drag(Figure 4.3 (b) and (c)). However in this particular case, a diminution of the lift-to-drag ratio corresponds to an increase of the lift coefficient and a faster increase of the drag coefficient.

Figure 4 .Figure 4

 44 Figure 4.4 shows the roll angle, the kite velocity, and the kite turn rate for the same trajectories and the same tether lengths as Figure 4.3. The results for the roll angles seem to be consistent with intuitive physics of kite: the large angles of roll are obtained in the case of the small radii of trajectory, and vice versa. For high ratio of radius over

  Figure 5.1 -Denomination of the main kite elements

Figure 5

 5 Figure 5.2 -Buckling of a LEI kite during flight. Picture taken during experimental tests of[START_REF] Behrel | Investigation of kites for auxiliary ship propulsion: experiment setup, trials, data analysis and kite specs novel identification approach[END_REF]

Figure 5

 5 Figure 5.3 -Real profile of a LEI kite. The black cross denotes the position of the sewing between the inflatable tube and the canopy. The blue curve represents the canopy in the IDK model.

Figure 5 . 4 -

 54 Figure 5.4 -Example of a LEI kite according to the IDK model

Figure 5 . 5 -

 55 Figure 5.5 -Kite as a Beam concept and division of the kite in several elementary cells

Figure 5 . 7 -

 57 Figure 5.7 -Geometry of the elementary cell

  Figure 5.8 -Shape of the canopy under homogeneous pressure loading with contour value representing the out-of-plane displacement along -→ z c (in meter)

Figure 5

 5 Figure 5.9 -Case (β): out-of-plane shear along z c . The contour value represents the displacement along z c .

Figure 5 .

 5 Figure 5.10 -Micro-beams and equivalent beams with the kite respective geometry

Figure 5 .

 5 Figure 5.11 -Determination of the centroid of a micro-beam

Figure 5 .

 5 Figure 5.13 -Convergence study of the D/CS ratio

  Figure 6.1 -Boundary conditions in displacements (U ) and angular velocity (VR) of the kite under pressure load

Figure 6

 6 Figure 6.3 -Description of the different comparison cases

Figure 6

 6 Figure 6.4 -Complex Finite Element model of the LEI kite.

  Figure 6.5 -Reaction forces and moment of the kite wing for the different comparison cases.

Figure 6

 6 Figure 6.6 -Complex Finite Element model of the cylindrical wing

  Figure 6.7 -Reaction forces and moment of the cylindrical wing for the different comparison cases.

Figure 6

 6 Figure 6.8 -Deformation of one batten during the first pressure step (in red) and after the displacement step (in green). Case at 200P a and displacement along the X axis.

Figure 6 . 9 -Figure 6 Figure 6

 6966 Figure 6.9 -Initial profile of the canopy

Figure 6 Figure 6

 66 Figure 6.12 -Evolution of the equivalent beam properties for position of maximum camber R LA ∈ [0.25; 0.5] for a range of pressure [100; 500P a], C MX = 0.065

Figure 6 . 14 -

 614 Figure 6.14 -Equivalent beam frame and micro-beam frame from the Kite as a Beam model of a kite geometry.

Figure 7 . 1 -

 71 Figure 7.1 -Set of bridles, tethers and connection to the wing for the IDK model (left) and the KaB model (right).

Figure 7 . 2 -SFigure 7

 727 Figure 7.2 -Description of the effective angle of incidence

  Figure 7.4 -Pressure coefficient provided by the XFOIL simulation on a LEI kite section at 4°of incidence

Figure 7 . 5 -Figure 7 Figure 7 . 8 -

 75778 Figure 7.5 -Boundary conditions on the IDK model for the structure iteration, global view (left) and close up (right).

Figure 7

 7 Figure 7.9 -Local twist of the aerodynamic sections after convergence for the wing at 7°ofincidence. The initial geometry twist is in gray.

Figure

  Figure 7.10 -Local effective angle of incidence and local projected wind velocity after convergence for the wing at 7°of incidence

Figure 7

 7 Figure 7.12 -Attachment point on the aerodynamic section of the KaB model (L KaB ) and the IDK model (L IDK )

Figure 7

 7 Figure 7.13 -Local twist of the kite wing for both model at 7°of incidence, after the first structure iteration in black and after convergence in red. (initial geometry twist in gray)

Figure 7 . 14 -

 714 Figure 7.14 -Wing at 7°on incidence after convergence with 20, 60, 100 or 150 sections over the span. (a) Magnitude of the aerodynamic loading as a function of the number of sections on the wing in the range [20; 150] for the FSI converged result in red and for the first iteration (with an undeformed wing) in black. (b) Relative difference on the aerodynamic loading with respect to the number of sections on the full span of the wing for the converged result in red and the first iteration in black.

Figure 7 .Figure 7

 77 Figure7.16 shows the deformation of the wing at 7°of incidence after convergence. The deformations of the wing only differ at the tips. The wing with the high mesh density closes more but has also less wingsweep.

Figure 7 .

 7 Figure 7.16 -Geometries of the wing at 7°of incidence after convergence with 20, 60 or 150 sections over the span. (a) Solid line represents the lifting line coordinates. (b) Solid line represents the leading edge, dashed line the trailing edge.

  (b)), but on the deformations of the wing. To take into account the influence of the work pressure, the equivalent beam properties could be computed for different work pressures, 3 for example, as Section 6.3 showed that the evolution of the equivalent beam properties are almost linear. The work pressure of each elementary cell can be estimated from the lifting line model and as a first approximation, a linear interpolation would give the equivalent beam properties of each elementary cell as a function of an updated work pressure

  Figure 7.17 -Wing at 7°of incidence after convergence with the IDK model (in red), the KaB model with a work pressure of 150P a (in blue), and with work pressures extracted from the IDK model (in green). (a) Local twist of the aerodynamic sections. (b) Magnitude of the local aerodynamic loading per unit length. (c) Solid line represent the leading edge, dashed line the trailing edge. (d) Deformed lifting line.

Figure 7 .

 7 Figure 7.18 -Wing at 7°of incidence after convergence with the IDK model (in red), the KaB model with a work pressure of 150P a (in blue), with calibrated properties (in green). (a) Local twist of the aerodynamic sections. (b) Magnitude of the local aerodynamic loading per unit length. (c) Solid line represent the leading edge, dashed line the trailing edge. (d) Deformed lifting line.

  The wingsweep is equal to d 1 at the tips of the wing with an elliptic distribution, but one can choose the slope at the tips of the tangent at the wingsweep curve f (s).

Figure A. 1 -

 1 Figure A.1 -Description of the lift polar curve

  ). . . . . . . . . . . . . . . . . . . . . . . 94 7.2 Description of the effective angle of incidence . . . . . . . . . . . . . . 96 7.3 Real section of a LEI kite (in red and blue) and XFOIL section (in blue, red and black). S J denotes the beginning of the added solid part in XFOIL, S S denotes the sewing point between the canopy and the inflatable tube. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 7.4 Pressure coefficient provided by the XFOIL simulation on a LEI kite section at 4°of incidence . . . . . . . . . . . . . . . . . . . . . . . . . . 97 7.5 Boundary conditions on the IDK model for the structure iteration, global view (left) and close up (right). . . . . . . . . . . . . . . . . . . . 99 7.6 Magnitude of the aerodynamic loading on the kite for each fluid-structure iteration step. The IDK FSI model was run for the case of a 7°of incidence with a wind velocity of 20.45m.s -1 , with another wind velocity of 12m.s -1 in blue and the KaB FSI model was tested for the case of 3°of incidence with a wind velocity of 20.45m.s -1 ( in black). . . . . . 99 7.7 Global deformations of the wing at 7°of incidence for the IDK model (in red) and for the KaB model (in blue). In gray the initial geometry. Solid lines represent the leading edges, dashed lines the trailing edges. 101 7.8 Deformations of the lifting line of the wing at Wing at 7°on incidence after convergence with 20, 60, 100 or 150 sections over the span.

7°of incidence for the IDK model (in red) and for the KaB model (in blue). In gray the initial geometry. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 7.9 Local twist of the aerodynamic sections after convergence for the wing at 7°of incidence. The initial geometry twist is in gray. . . . . . . . . . 102 7.10 Local effective angle of incidence and local projected wind velocity after convergence for the wing at 7°of incidence . . . . . . . . . . . . 103 7.11 Magnitude of the local aerodynamic loading per unit length after convergence for the wing at 7°of incidence. . . . . . . . . . . . . . . . . . 103 7.12 Attachment point on the aerodynamic section of the KaB model (L KaB ) and the IDK model (L IDK ) . . . . . . . . . . . . . . . . . . . . . . . . . 104 7.13 Local twist of the kite wing for both model at 7°of incidence, after the first structure iteration in black and after convergence in red. (initial geometry twist in gray) . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 7.14

  Converged FSI results of the kite wing at 7°of incidence for the IDK model, for the KaB model with a work pressure of 150P a (20 Sections) and for the KaB model with a work pressure in the range [250; 300P a] extracted from the IDK results (IDK pressure). . . . . . . . . . . . . . . 7.3 Converged FSI results of the kite wing at 7°of incidence for the IDK model, for the KaB model with a work pressure of 150P a (KaB 20 Sections) and for the KaB model with calibrated properties (KaB Calibrated). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

between the IDK and the KaB models for the cylindrical wing. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.1 Global results of the two FSI simulations at 7°of incidence . . . . . . 7.2 xvii Symbols used

Table 2 .

 2 1 -Number of iterations for each of the initial circulations.

Table 2

 2

	.2 and will be

Table 2 .

 2 2 -Estimation of the confidence interval for a number of sections n between [30; 100], α is the angle of incidence of the wing

  It can also be noticed inFigures 3.7 and 3.8 (c) that intermediate meshes in the buffer layer, for y+ in the range of[5; 30], led as expected to highly model dependent results, and that the Spalart-Allmaras model is valid only for fine meshes around y+ = 1.

					Figure (c)
		Figure (a) 2Crms/Ccvg 2Crms/Ccvg 2Crms/Ccvg 2Crms/Ccvg Figure (b) y+ ∼ 1 y+ ∼ 30	Cdif
	Lift	0.6%	0.7%	0.8%	1.2%	1.2%
	Drag	1.7%	2.5%	5.0%	1.5%	2.0%
	Moment	0.9%	1.5%	1.5%	2.4%	2.3%

Table 3 .

 3 1 -Estimation of the confidence intervals, coming from the variations of the domain size, of the mesh and of the turbulence model. C stands for coefficient.

Table 3 .

 3 .2. 2 -Overview of the different validation cases

	Chord	Sweep angle Twist angle
	case 1 linear	null	null
	case 2 linear	constant (30°) null
	case 3 non-linear non-linear	Linear

Table 6 .

 6 1 -Boundary conditions in displacements (U ), velocity (V ) and angular velocity (VR) on the IDK model for the different comparison cases

Table 6 .

 6 2 -Boundary conditions in displacements (U ), rotation (UR) and angular velocity (VR) on the KaB model for the different comparison cases

Table 7 .

 7 2 -Converged FSI results of the kite wing at 7°of incidence for the IDK model, for the KaB model with a work pressure of 150P a (20 Sections) and for the KaB model with a work pressure in the range [250; 300P a] extracted from the IDK results (IDK pressure).

	929	390	400

Table 7 .

 7 KaB 20 Sections KaB Calibrated Final position (°) 3 -Converged FSI results of the kite wing at 7°of incidence for the IDK model, for the KaB model with a work pressure of 150P a (KaB 20 Sections) and for the KaB model with calibrated properties (KaB Calibrated).

		8	12	9
	Aerodynamic loading (N)	929	390	580
	Lift coefficient	0.712	0.271	0.432
	Drag coefficient	0.378	0.130	0.121
	Wingsweep displacement (m) Opening displacement (m) Front tether loading (N)	-0.249 0.13 560	-0.407 -0.02 187	-0.344 0.52 249
	Back tether loading (N)	365	202	333

U Y (P ) = (P -C). #» X tan(ω) (6.4)where ω is the angle of rotation along Z axis. The displacement U Z is set to zero at the

Remerciements

The Young modulus is given by:

(5.3)

Which leads to the correct bending stiffness:

The shear modulus is set as:

(5.5) and the torsion stiffness as:

The transverse shear stiffness of the section is computed according to Cowper [START_REF] Cowper | The shear coefficient in Timoshenko's beam theory[END_REF]) as:

(5.7)

Abaqus™ defines a Slenderness Compensation Factor (SCF) (SIMULIA (2014)) to prevent the shear stiffness from becoming too large in slender beam element. The transverse shear stiffness of a linear beam element of length L is defined in Equation 5.8:

(5.8)

In order to obtain the transverse shear stiffness of a Timoshenko beam of length L, SCF has been set to:

The option *GENERAL SECTION of Abaqus™ allows also to split the centroid of the beam section from the element which can be useful in some instances (see Section 5.3). model computation time (see table 0.3). Actually, there is more than a factor 100 in the computation time between the two models , keeping in mind that the inflation step is included in the IDK simulation time. If the computation of the equivalent beams properties is taken into account, the comparison is less favorable for the KaB model, with a computation time of about 40s per cell. However, these computations are conducted just once for a given pressure. The used range of pressure (from 100P a to 400P a) corresponds to classic flight conditions [START_REF] Behrel | Investigation of kites for auxiliary ship propulsion: experiment setup, trials, data analysis and kite specs novel identification approach[END_REF]) and is in the intended use range, considering a maximal wind loading for the use of kites of 100kg/m 2 .

Comparison for a LEI kite geometry

Firstly the two models are compared on a classic LEI kite geometry. The chosen kite was designed by beyond the sea® for a surface of 50m 2 , see Figure 6.4. The geometry is the same for the two models beside the geometrical approximation made to build the elementary cell of the KaB model. The kite has a mean radius of 4.3m, the chord length varies from 4m at the root to 1.3m at the tips with a maximal sweep of 0.36m.

The radius of the leading edge and the inflatable battens vary in the same way for the two models, the material properties of the canopy and the battens are identical and modeled with the same finite elements. The elastic behavior of the canopy fabrics is isotropic while the behavior of the inflatable tubes is orthotropic.

Simulation results can be seen in Figure 6.5 and the difference between the two models computed as described in Equation 6.6 can be seen in Table 6. Table 6.4 -Gap (in %) between the IDK and the KaB models for the kite geometry displacement.

Global results

The results presented here are those after convergence for the wing at 7°of incidence with a wind of 20.45m/s . The global results are presented in Table 7.1. The final position is the angle in the plane (O; #» x ; #» z ) between #» z and the front tether. The wingsweep is the displacement along the #» x axis at the tips. The opening is the displacement of the tips of the wing in #» y direction, that means the distortions of the wing along the #» y axis. The tether loading is the magnitude of the reaction force vector at the attachment points on the ground. The relative difference between the two models is given by taking the IDK FSI as reference.

IDK model KaB model Final position (°)

8 34

Aerodynamic loading (N) Front tether loading (N)

Back tether loading (N)

365 91

Table 7.1 -Global results of the two FSI simulations at 7°of incidence

The differences between the two models are really significant. The relative difference for the aerodynamic loading is 66%. The initial loading for the undeformed wing is 560N , thus according to the IDK model, the wing is distorted and hence increases the aerodynamic loading, while it is the contrary for the KaB model. The same order of magnitude can be found for the difference in position in the plane (O; #» x ; #» z ), for the lift and drag coefficient or for the loading of the tethers, which is linked to the position of the kite in the symmetry plane (O; #» x ; #» z ). The geometrical deformations can also be seen in Figures 7.7 and 7.8. Both models predict the same direction in wingsweep displacement, but the KaB estimation is way larger (by 50%). The opening of the kite, the distance between the tips of the wing, is equal to 2.8m non-deformed which represents an evolution of 4.6% for the IDK model. If the wing opens in the IDK model, it closes in the KaB model, of 5.7%. Therefore it seems that the models agree only for the direction of the wingsweep displacement.

The global deformations of the wing can be seen in Figures 7.7 and 7.8. The leading edge (solid lines) and the trailing edge (dashed lines) are represented here, the IDK model in red, the KaB model in blue and the initial geometry in gray. In Figure 7.8, only the lifting line is represented for the sake of clarity. The wingsweep can be seen again in Figure 7.7 as well as the opening of the wing in Figure 7.8. The variations of the aerodynamic sections length is visible in Figure 7.7 with the deformations of the

Appendix

Appendix: Lifting Line model

Chord law

Different chord laws c(s) are implemented in the program. c r is the chord length at the root, c t the chord length at the tips. s is the curvilinear abscissa over the span and s max is the maximal curvilinear abscissa.

• Linear:

• Elliptic:

Twist law

Two twist laws α(s) are possible in the program. α r is the twist at the root of the wing, α t at the tips. s max is again the maximal curvilinear abscissa.

• Linear:

• 3rd degree polynomial: The twist is equal to α r at the root and α t at the tips of the wing. But one can set the twist at c v s max with the value

Wingsweep law

Three wingsweep laws f (s) are implemented in the program. d 1 is the wingsweep at the tip of the wing. s is the curvilinear abscissa over the span and s max is the maximal curvilinear abscissa.

Drag coefficient

The drag polar curve (see Figure A.2) is divided into two parts, part I is polynomial of degree 3 and part II is linear. The limit between them is given by α lim . The polynomial go through the points (α i , C di ) and (α min , C dmin ) is the minimum value of this part of the drag polar curve.

-10 0 10 20 0.00 0.50

I II

Angle of incidence (°) Drag Coefficient

Moment coefficient

The moment polar curve (see Figure A.3) is also divided into two parts, but both polynomial of degree 3 for part I and of degree 2 for part II. α lim is the limit between the two parts of the curve. The first part pass through the points (α i , C mi ) with i ∈ [1; 4] and the second part through the points (α i , C mi ) with i ∈ [5; 7].

-10 0 10 20 -0.05 0.00 0.05

I II • For the canopy: Young modulus (N/m 2 ), Poisson's ratio and grammage of the fabric (g/m 2 )

Angle of incidence (°) Moment Coefficient

• E: Fictive thickness corresponding to the following mechanical properties of the fabrics (m)

• For the leading edge tube : Longitudinal modulus E mc and shear modulus G mc in N/m. Radius of the leading edge for each elementary cell (m)

• For the battens: Longitudinal modulus E mc and shear modulus G mc in N/m. Minimum and maximum radius for each elementary cell (m)

• For the leech line: Density (kg/m 3 ), diameter (m) and Young modulus (N/m 2 )

• Attachment points of the bridles over the span of the kite

• Diameter of the front tethers (m)

• Diameter of the back tethers (m)

• Density of the tethers (kg/m 3 ) and Young modulus (N/m 2 )

• Position on the ground of the back tether (m). Front tethers are set in (0; 0; 0).

• Lengths of the various bridles (m)

• Lengths of the front and back tethers (m) Titre : Modélisation avec prise en compte de l'interaction fluide-structure du comportement sous charge d'un cerf-volant pour la traction auxiliaire des navires Abstract : The present thesis is part of the beyond the sea project which aims to develop tethered kite systems as auxiliary devices for ship propulsion. As a kite is a flexible structure, fluid-structure interaction has to be taken into account to calculate the flying shape and aerodynamic performances of the wing. A 3D Non-Linear Lifting Line model has been developed to deal with non-straight kite wings, with dihedral and sweep angles variable along the span and take into account the non-linearity of the section lift coefficient. The model has been checked with 3D RANSE simulations over various geometries and produces satisfactory results for range of incidence and sideslip up to 15°, with typical relative differences of few percent for the overall lift. The local results are also correctly estimated, the model is able to predict the position of the minimum and maximum loading along the span, even for a wing in sideslip.

Simultaneously, a structure model has been developed. The core idea of the Kite as a Beam model is to approximate a Leading Edge Inflatable kite by an assembly of beam elements, equivalent each to a part of the kite composed of a portion of the inflatable leading edge, two inflatable battens and the corresponding canopy.

The Kite as a Beam model has been compared to a complete kite Finite Element model over elementary comparison cases. The results show the behaviour differences of the two models, for example the torsion stiffness is globally overestimated by the Kite as a Beam model. Eventually, the Kite as a Beam model coupled with the 3D Non-Linear Lifting Line model is compared to the complete Finite Element model coupled with the 3D Non-Linear Lifting Line model. The gain in computation time is really significant but the results show the necessity of model calibration if the Kite as a Beam model should be used to predict the results of the complete Finite Element model.