Je souhaite remercier mes directeurs de thèse Angelo Iollo et Michel Bergmann. J'ai beaucoup apprecié de travailler avec vous pendant ces années de thèse. Nos échanges professionnels et personnels ont été très enrichissants.

Je remercie vivement Kai Schneider et Éric Lamballais pour avoir accepté de rapporter mon mémoire de thèse. Vos remarques ont été très précieuses pour améliorer le manuscrit de la thèse. Je remercie également Charles-Henri Bruneau pour avoir présidé ce jury de thèse, ainsi que Haysam Telib, Jean-Yves Grandidier et Andrea Ferrero

Contents

Introduction Chapter 1.

Wind Turbine Aerodynamics and Existing Modelling Tools

Introduction Industrial Context

Almost all states worldwide agreed during the 2015 United Nations Climate Change Conference that was held in Paris on the reduction of climate change. One of the goal is the limiting of global warming to "well below 2°C" compared to pre-industrial levels. The agreement calls also for zero net anthropogenic greenhouse gas emissions to be reached during the second half of the 21 st century. This conference didn't envisage concrete actions to implement, no sanction is moreover planned in the case where a country doesn't make any eort. However, this agreement proves that the world is realizing the importance to limit climate change. To reach this objective of reducing greenhouse gas emissions, the development of renewable energies will play a fundamental role and a likely winning scenario is a mix of several sources of energy (wind, solar, hydraulic, methanation, ...). Inside this mix, wind energy has several signicant advantages that will certainly lead to its further development in the next years. Wind is a local resource and its cost is becoming cheaper and cheaper thanks to technological improvements. Moreover, unlike the solar energy, wind can be available in a more regular way all along the day and night. France has also the advantage of having several main windy corridors that allow to have almost always somewhere the wind resource available. The possible extension to oshore wind plants could also lead to a signicant development of wind energy. Indeed, the impacts are lower than onshore, while having a constant and regular wind. This could occur in the near future since costs of installations are constantly decreasing and the cost of oshore wind is currently lower than the new generation of nuclear power plant (EPR) under development in Flamanville (France).

This development of the wind energy sector is thus leading to technological advances. It can be seen in gure 1 that shows the evolution of the size of wind turbines over time. Thanks to the large blades, wind turbines are now able to catch wind even in areas where the wind resource is quite small. Concerning Valorem, the company is involved in the whole value chain of renewable energy with the development part, the construction of plants, the operation of wind farms and maintenance. The development of a wind plant consists in identifying the potential areas, the associated constraints and in obtaining the administrative authorization. Also, at this stage, some studies should be done to have an estimation of the wind resources and of the future energy production (as a function of layout of wind turbines and wind turbines themselves).

Environmental studies should also be done to estimate the impact of the future wind farm on the landscape and the fauna. Finally, electrical connection studies need to be integrated into the development phase to plan the grid connection of the wind farm. Valorem owns also several power plants (solar, wind or hydroelectric). Concerning the wind sector, which is the major energy developed at Valorem, it is essential to predict in an accurate way the future production of the wind plant. Indeed, the business model of the future wind farm is highly dependent on that element and a valid business model is necessary for the viability of Valorem. Therefore, the use of adequate tools is essential for the wind assessment studies. For a company as Valorem, The simulation of the eects of a gust on a wind turbine is a fundamental step in the determination of the maximum possible load encountered by the wind turbine during its life cycle. Indeed, the magnitude of a gust can be suciently high to increase the loads up the point of fracture. Catastrophic failures are reported in the literature for both wind turbine blades and towers (see for example [START_REF] Chou | Failure analysis and risk management of a collapsed large wind turbine tower[END_REF] and [START_REF] Chou | Failure analysis of wind turbine blade under critical wind loads[END_REF]). Therefore, Valorem was part of the European Horizon 2020 project AeroGust ("Aeroelastic Gust Modelling"). The aim of the project is to better understand the eects of wind gusts on both the aircraft and on wind turbine blades. As lifetime extension of the blades is of high interest for Valorem, being part of that project was a great opportunity.

Background to the Research

Since the Reynolds number of the involved bodies is quite large (with an order of several millions) and the ow eld is fully turbulent, a large range of turbulent structures is present. In order to achieve a good representation of the ow by using a reasonable amount of computational power, a turbulence model is required. Recently, several promising simulations of the ow eld around wind turbines have been carried out by using Large Eddy Simulations (see for example [START_REF] Arjomandi | Large eddy simulation of the wind turbine wake characteristics in the numerical wind tunnel model[END_REF], [START_REF] Norris | Wind turbine wake modelling using Large Eddy Simulation[END_REF] and [START_REF] Lorenzo | Validation of four LES and a vortex model against stereo-PIV measurements in the near wake of an actuator disc and a wind turbine[END_REF]). This is the path which will be followed in the present work devoted to the simulation of a full wind turbine. Finally, very few existing numerical tools are able to compute these analyses. Moreover, a key parameter is the time needed by engineers to dene an appropriate grid discretization for a given problem. This is a major issue explaining why INRIA and VALOREM chose to develop an immersed boundary code based on octree discretization where the denition of a mesh is straightforward.

coupling with the structural model of the blade.

The aim of the present thesis is to develop a rst base of a numerical tool that will allow to better understand the behaviour of wind blades. This tool could then be improved and adapted to dierent applications as wished by Valorem.

Outline of the Thesis

The rst chapter describes the basic aerodynamics of wind turbine to highlight the data of interest to compute. The existing numerical models used for preliminary design are then discussed.

A particular attention is brought on the Blade Element Momentum theory which is widely used in industrial process. Some numerical models in Computational Fluid Dynamics are then considered. We focused on those used for wind turbine modelling as it is our nal application, and on existing schemes that are based on an octree discretization. Dierent existing computational grids have indeed been studied and the octree conguration shows several advantages for our application. Advection schemes have also been studied since it led to an issue as realized during the thesis.

All the numerical methods developed during this work are explained in the chapter 2. Firstly, the choice of an immersed boundary method is commentated and the solution of Navier-Stokes equations in that context is described. The two-dimensional nite volume solver developed with quadtree discretization is then widely explained with the chosen stencils and interpolation methods. Each operator implemented in the code in order to solve the Navier-Stokes equations is validated thanks to the use of analytical functions and the order of convergence is computed. The extension of these original methods to three dimensional case follows with the same detailed explanations and the operators are preliminary validated by the use of analytical functions and the computation of the order of convergence. To end this chapter, the numerical models used to deal with high-Reynolds ow in our work are presented. Large Eddy Simulation with a Vreman subgrid model has been implemented to model turbulence. Wall functions are also employed to model the behaviour of the ow around the obstacles. Indeed, the octree discretization prevents us to be able to properly model the boundary layer near walls.

little experimental data is indeed available with the geometry of blades, the structural model and measurements on the blade. And none of this correlates that with the incoming wind. The context of the European project AeroGust1 allowed to collect this experimental data on an in-service wind turbine. All the specications needed to conduct the experimental work have been determined. Dierent technical solutions have been investigated to measure both the wind and to get data from the blade and the nal choice is presented. Then, on-site installation of the sensors has been planed and organized. Finally, the experimental data collected has been studied in order to better understand the eects of wind on the blade.

The nal chapter describes the application of the solver developed which is the simulation of the ow around a wind turbine blade. Since we have available experimental deformation data of the wind blade APX48, the structural model of this blade has been also implemented in the octree solver. The way the coupling between the uid and the structure is done in present work is also described. Studying the ow around a whole wind turbine required too many computational time. Therefore, to reduce that, it was chosen in the thesis to study only one blade and moreover, to use a rotating framework to focus the computational domain on the wind blade. An interesting feature of our developed solver on octree is the possibility to perform Adaptive Mesh Renement of the grid to optimize the grid size while being as precise as possible on the wake of the rotating wind blade. The implementation of this process is described together with the preliminary results obtained. Finally, all the features implemented have been put together to perform simulations of the ow over a rotating blade with real operating conditions of the wind turbine.

Chapter 1

Wind Turbine Aerodynamics and Existing Modelling Tools

This chapter explains the operation of a wind turbine with the aerodynamic forces that act on the blades. In order to compute these quantities and estimate the power extracted by some wind turbine, several kind of modelling exist. Simplied numerical models can be used for preliminary design. These models are mostly based on the Blade Element Momentum theory as described in [START_REF] Mikkelsen | Actuator Disk Methods Applied to Wind Turbines[END_REF][START_REF] Sorensen | Unsteady actuator disc model for horizontal axis wind turbine[END_REF][START_REF] Sanderse | Aerodynamics of wind turbine wakes -Literature review[END_REF]. These models are bi-dimensional whereas the ow around wind turbines is intrinsically tri-dimensional. Therefore, a signicant work of modelling is necessary to accurately reproduce the physics of the ow with such models. An other approach is the Computational Fluid Dynamics (CFD) that rest on the Navier-Stokes equations and their approximation. As this subject is very wide, a bibliographic study has been performed in order to highlight the numerical schemes adapted to the problem of the ow around a rotating wind turbine.

Wind Turbine Aerodynamics

Aerodynamics of the blades has a signicant importance on the energy eciency of the wind turbines. The aerodynamic forces exerted on the blades can be seen in gure 1.1. U is the incoming wind speed and U ω corresponds to velocity resulting from the rotational speed of the wind turbine. So, composition of wind speed and rotational speed is named relative speed (U r ) and corresponds to the velocity arriving on the blade. At blade root, the rotational speed is low which implies that the relative speed is really close to the wind speed. On the contrary, at tip blade, the rotational speed is signicantly higher than the wind.

The parameter β represents the sum of the twist angle at considered radius and of the pitch angle, which is the orientation angle of the blade around its axis (this angle varies along time). β, named the angle of incidence is the angle between the chord line and the rotational plane of the rotor. α is the angle of attack, being the angle between the chord line and the relative speed (U r ). The ow angle, Φ can then be calculated as: Φ = α + β.

Due to the air ow around the airfoil, two forces exerted on the blade: drag (D) that act opposite to the relative ow and lift (L). Lift is the component perpendicular to the oncoming ow direction and corresponding to the depression due to acceleration of the uid on the suction side of the airfoil. The sum of these two vectors forms the aerodynamic force F and multiplying it with the distance to the rotor's center point allows to get the torque. Its tangential component (F T ) contributes to the rotation of the wind turbine and the normal component (or axial component, referred to the rotor axis) of the aerodynamic force (F N ) corresponds to the thrust,

U U r -U ω L D F N F T α β Φ Figure 1
.1: Sketch of the aerodynamic forces exerted by the ow around a wind blade airfoil that subjects the rotor to mechanical stresses. The computation of these aerodynamic forces is necessary to evaluate the energy eciency of the wind turbine.

The power extracted by a wind turbine is a signicant variable to analyse and the power coefcient can be dened as the ratio between power extracted by the turbine and available wind power:

C P = P extract P available = P 1 2 ρAV 3 , (1.1)
with A = πD 2 4 , the rotor swept area. In the same way, the thrust coecient (notion dened above) can be introduced:

C T h = F N 1 2 ρAV 2 .
(1.

2)

The power coecient C P depends on the rotational speed of the rotor and can be expressed as a function of the tip speed ratio λ that represents the ratio between rotational speed at tip blade and the wind speed:

λ = ΩR U , (1.3)
Ω being the rotational speed of the rotor (in rad.s -1 ) and R the blade radius. The wind turbine torque T , as previously seen can be written:

T = P Ω = r ∧ F , (1.4)
with r being the position vector (vector from the origin of the coordinate system dened, to the point where the force is applied) and F is the aerodynamic force vector. The C Pλ curves are used in wind turbine design to determine the rotor power for any combination of wind and rotor speed. In gure 1.2 can be seen the corresponding curve for the Euros EU120 wind blade. For this blade, the maximal power coecient is around 0.49 and is obtained for λ opt = 9.1. That means, so, for a blade tip speed equal to 9.1 times wind speed. Using this curve as well as data from the wind turbine (its rated capacity, its minimal and maximal wind speed, ...) it is possible to retrieve the power curve of the wind turbine. The power curve of a wind turbine is a graph that indicates how large the electrical power output will be at dierent wind speeds. An example can be seen in gure 1. 3.

Introducing the well-known Mach number as the ratio of ow velocity to the speed of sound, the incompressibility condition of ow around a wind turbine can be studied. Indeed, compressibility eects are small if the Mach number is less than 0.3. The speed of sound is equal approximately to 343 m.s -1 for air at 20°C. The maximal wind speed around a wind turbine appears at tip blade. Depending on the wind turbine model, its value is ranged between 70 m.s -1 for some 2 M W wind turbine and until 100 m.s -1 for some oshore 5 M W wind turbine. That implies Mach number included in the range 0.2 -0.3 for nominal tip speeds which supports the incompressibility assumption for the simulation of the ow around a wind turbine.

Existing Simplied Numerical Models for Preliminary Design

The existing tools to simulate wind turbine aerodynamics can be divided in two main categories: Computational Fluid Dynamics (CFD) and simplied models (or hybrid models). CFD computes the whole Navier-Stokes equations in a three-dimensional domain containing the wind turbine and evaluates the aerodynamic forces on the blades dynamically. Hybrid models replace the real wind turbine geometry by simplied models of the ow around the wind turbines. All the hybrid models are based on blade element momentum theory to compute the forces on the wind blades and are coupled with Reynolds' averaged Navier-Stokes equations (RANS). This section explains the blade element momentum (BEM) theory before describing actuator disc and actuator line models. These two models are named the way the blade geometry is replaced with.

Blade Element Momentum (BEM)

BEM theory is one of the oldest aerodynamic model in the eld of wind engenering and is widely used in industrial world for the design of wind blades [START_REF] Sorensen | Aerodynamic aspects of wind energy conversion[END_REF][START_REF] Moriarty | AeroDyn theory manual[END_REF][START_REF] Bot | Blade Optimization Tool (BOT): User Manual[END_REF]. This model is a combination of two theories: axial momentum and blade element theory. This method has been developed by Betz in 1926 and is explained more recently in many articles [START_REF] Mikkelsen | Actuator Disk Methods Applied to Wind Turbines[END_REF][START_REF] Sorensen | Unsteady actuator disc model for horizontal axis wind turbine[END_REF][START_REF] Sanderse | Aerodynamics of wind turbine wakes -Literature review[END_REF]. Basically, in axial momentum theory, a unidimensional model is used to express the wind extracted energy and the ow is supposed to be incompressible. The blade element theory takes into account the uid's rotational component that comes from blade's rotation. These theories are detailed below.

Axial momentum theory

In this theory, only an axial, incompressible ow is assumed, with airow passing through the rotor without friction. The modelling consists in considering the wind turbine's rotor as a permeable disc without thickness which is set in motion by the airow passing through it. As shown in gure 1.4, the air passing through the wind turbine rotor is subjected to a change of speed. The axial induction factor (often used in literature) is introduced and corresponds to the wind speed decrease fraction between upstream value V 1 and the value of wind speed going through rotor's plane V 2 :

a = 1 - V 2 V 1 .
(1.5)

The variable V 4 corresponds to the downstream wind speed. Using the law of conservation of momentum, the theory considers that the force acting on the actuator disc reduces to the variation of axial momentum between upstream and downstream sections of the stream tube:

F N = ṁ(V 1 -V 4 ). (1.6)
Here, the mass ow through the rotor can be written ṁ = ρAV 2 , with A being the rotor area. In the same manner, using conservation of energy, the energy extracted by the disc can be written as:

E = 1 2 m(V 2 1 -V 2 4
).

(1.7)

Figure 1.4: Axial stream tube around a wind turbine [START_REF] Ingram | Wind turbine blade analysis using the Blade Element Momentum method[END_REF] The power extracted from the oweld by the rotor can thus be written in two manners. On the one hand, the power, by denition, can be written as the energy transferred per unit time. Using eq. (1.7), we nally have:

P = 1 2 ṁ(V 2 1 -V 2 4
).

(1.8)

On the other hand, from the aerodynamics, with the axial ow condition, we have the denition:

P = F N V 2 .
With the use of eq. (1.6), we get:

P = ṁ(V 1 -V 4 )V 2 .
(1.9)

From eqs (1.9) and (1.8), can be deduced:

V 2 = 1 2 (V 1 + V 4 ).
(1.10)

By replacing V 4 in eq. (1.6), we get:

F N = 2ρAV 2 1 a(1 -a).
(1.11)

Energy eciency

With eqs (1.1), (1.5), (1.9) and (1.10), the non-dimensional power coecient is established as:

C P = ρAV 2 2 (V 1 -V 4 ) 1 2 ρAV 3 2 = 4a(1 -a) 2 .
(1.12)

By solving dC P da = 0 for eq. (1.12), the optimal conversion of energy possible can be found for a = 1 3 , which corresponds to C Pmax = 16 27 0.59. The highest energy eciency of a ideal wind turbine is so around 59%. This value is usually reered to as the Betz limit. The axial momentum theory is a simple unidimensional model that doesn't take into account the air rotational ow.

Blade element theory

Blade element theory assumes an axial upstream ow far from rotor but adds the rotational speed of the ow at rotor's plane level. This model considers that the blade is constituted of several elements. The theory relies on the assumption that there is no aerodynamic interaction between dierent elements, so each element is independent from the others. Figure 1.5 shows an Figure 1.5: Blade element theory [START_REF] Ingram | Wind turbine blade analysis using the Blade Element Momentum method[END_REF].

element with a thickness of dr that constitutes a blade and ctious ring formed during blade's rotation. An other assumption of the theory is that the ow around each element is only bidirectional, thus, interaction forces on a blade element reduces to lift and drag. After computing the forces on each element, the global blade performance is evaluated by integrating the stresses exerted on each element.

The tangential ow induction factor a is introduced in blade element theory. Its denition is a = ω Ω , with ω the angular velocity of the ow on the rotor's plane and Ω, the rotational speed of the rotor. Using notation of gures 1.1 and 1.5, these both relations can be written, according to the denition of a (see above) and a :

U = U ∞ (1 -a),
(1.13)

U ω = Ω + ωr = Ωr(1 + a ), (1.14)
with U ∞ that corresponds to the far upstream wind speed. Then, geometrically, the following relations can be written:

tan(Φ) = U ∞(1-a) Ωr(1 + a ) , (1.15)
U r = U 2 ∞ (1 -a) 2 + Ω 2 r 2 (1 + a ) 2 .
(1.16)

The basic hypothesis of the theory is that the force exerted on a blade element equals variation of momentum of the air passing through the blade. Thus, lift and drag are:

dL = 1 2 ρU 2 r BC L dS, (1.17) dD = 1 2 ρU 2 r BC D dS, (1.18)
with C L and C D being respectively the lift and drag coecients. dS = cdr is the surface of the considered blade element and B corresponds to the number of blades. By projecting, the tangential force (related to torque) acting on a blade element can be obtained:

dF T = 1 2 ρcdrU 2 r B(C L sin(Φ) -C D cos(Φ)). (1.19)
The torque exerted on a blade element (dT ) being the tangential force time the radius, we can write:

dT = 1 2 ρcrdrU 2 r B(C L sin(Φ) -C D cos(Φ)). (1.20)
This is the force that need to be computed in order to get the power developed by the wind turbine. In the same way, the normal force (related to thrust) can be written:

dF N = 1 2 ρcdrU 2 r B(C L cos(Φ) + C D sin(Φ)). (1.21)
By introducing the denition of local solidity σ = Bc 2πr , eqs (1.20) and (1.21) become:

dT = σπρr 2 drU 2 r (C L sin(Φ) -C D cos(Φ)), (1.22
)

dF N = σπρrdrU 2 r (C L cos(Φ) + C D sin(Φ)). (1.23)
Combination of both theories -Blade element momentum

Eq. (1.11) can be rewritten with notation from blade element and taking an area of the blade element A = 2πrdr allows to obtain:

dF N = 4ρπrdrU 2 ∞ a(1 -a).
(1.24)

In the same way, using notation of gure 1.4 and extending axial momentum theory, variation of momentum in tangential direction gives us the tangential force that is exerted on the blade:

dF T = ρAV 2 V T , (1.25)
with V T = ω 3 r, V T being the tangential velocity of ow just after the rotor and ω 3 is dened as the angular velocity. In order to determine ω 3 , the conservation law of angular momentum is applied between section 1 and 2 of gure 1.4.

ω 1 r 2 1 = ω 2 r 2 2 .
(1.26)

As the upstream ow far from the turbine is only axial, we have ω 1 = 0. That implies from eq.

(1.26) that ω 2 = 0. It can be deduced that angular velocity is subjected to a discontinuity at rotor's plane, this is due to the torque applied on the rotor. The angular velocity of the ow on the rotor's plane is the average of those just before and after the rotor.

ω = ω 2 + ω 3 2 = ω 3 2 . (1.27)
Finally, using denition of a and a , eq. (1.25) becomes:

dF T = ρ(2πrdr)V 1 (1 -a)2a Ω.
(1.28)

We obtain then the torque, noting that V 1 corresponds to U ∞ using notation of gure 1.1.

dT = 4ρU ∞ (1 -a)a Ωr 3 πdr. (1.29)
With both theories, torque and thrust have been expressed in two dierent forms. We can thus solve the problem and obtain with eqs (1.23) and (1.24):

a = 1 4sin 2 (Φ) σ(C L cos(Φ) + C D sin(Φ)) + 1
.

(1.30)

Use of eqs (1.22) and (1.29) allows to get:

a = 1 4sin(Φ)cos(Φ) σ(C L sin(Φ) -C D cos(Φ)) + 1 . (1.31)
Blade element momentum theory is thus a simplied model that allows to compute the forces that are applied on a rotor and deduce the power extracted from the rotor. The data needed for BEM are the following:

radius of blade R,
rotational speed of the rotor Ω, number of blades B, distribution of airfoils with associate chord length c and twist angle β, polar curves of each airfoil along the blade, that allows to know drag and lift coecient of each airfoil.

As no analytical expression exists for the ow induction factors a and a , an iterative process is used to determine them. Convergence of results of a and a means that axial momentum and blade element theories converge.

Generalized Actuator Disc Model

Generalized actuator disc (GAD) method can be seen as an extension of BEM theory. Indeed, the wind turbine's rotor is modeled as a zero thickness disc in a uid ow. But, instead of considering each blade element independent from each other as in BEM, GAD model computes Navier-Stokes equations to get the velocity eld, especially on the blades. The aerodynamic forces are then obtained with BEM method as described above. Axisymmetric versions of Navier-Stokes equations are developed for GAD model. Numerical implementation can be based on nite difference as developed by Sorensen and Mikken [START_REF] Sorensen | Unsteady actuator disc model for horizontal axis wind turbine[END_REF] and presented by Mikkelsen [START_REF] Mikkelsen | Actuator Disk Methods Applied to Wind Turbines[END_REF]. The model consists in supposing an axisymmetric tube, thus, a 2D plane (r, z) is considered, discretized with Cartesian mesh, where Navier-Stokes equations are solved.

Previously to this thesis, a work has been done with development of two codes based respectively on GAD and BEM theory. These models have been compared with experimental results obtained in a wind tunnel on a two blade wind turbine of 5 m diameter. The experimental process is described in [START_REF] Hand | Unsteady aerodynamics experiment phase VI: Wind tunnel test congurations and available data campaigns[END_REF]. The predicted power as a function of wind speed computed with BEM and GAD models have been compared with these experimental results and can be seen in gure 1.6. The results show that for high wind speeds, both models are inaccurate. This explains why an other approach with full 3D Computational Fluid Dynamics (CFD) has been studied.

Actuator Lines

The actuator line model consists in replacing the actual geometry of the rotor blades by lines carrying body forces corresponding to the loading of rotor blades. This technique was originally developed by Sørensen and Shen [START_REF] Sorensen | Numerical modelling of wind turbine wakes[END_REF] and presented later in [START_REF] Sorensen | Simulation of wind turbine wakes using the actuator line technique[END_REF]. The method is based on the same ideas as GAD. Dierences are, rst, the solving of three-dimensional Navier-Stokes equations in the domain. Large Eddy Simulation is employed, that applies a low-pass lter on Figure 1.6: Comparison of predicted power for BEM and GAD models with experimental data of NREL UAE phase 6 (experimental process described in [START_REF] Hand | Unsteady aerodynamics experiment phase VI: Wind tunnel test congurations and available data campaigns[END_REF]).

the equations. So, the scales below the grid scale are modelled with a sub-grid model. Secondly, the inuence of the wind turbine is modelled with lines instead of elements. The actuator lines are fully coupled with the CFD solver so that they can deect and rotate. This allows an interesting simulation of the rotor's wake. However, as in GAD, the forces applying on the blade are obtained with the polar curves that gives the drag and lift coecients for a well-known airfoil geometry as a function of the angle of attack. These models are thus adapted only when simulating an existing blade with a well-known airfoil distribution.

Computational Grid

Dierent categories of grids exist to discretize Navier-Stokes equations. Structured grids are characterized by regular connectivity with quadrilateral elements in 2D and hexahedra in 3D. This kind of mesh has a better convergence and higher resolution compared to unstructured grids [START_REF] Castillo | Mathematical Aspects of Numerical Grid Generation[END_REF][START_REF] George | Automatic mesh generation: application to nite element methods[END_REF]. Unstructured meshes typically employ triangles in 2D and tetrahedra in 3D. They are well adapted for complex geometries so that the mesh can t the obstacles. Finally, hybrid grids exist that are mixing structured and unstructured blocks.

The octree data structure

The octree structure can be viewed in its dual nature, the tree and the grid, as it is represented in gure 1.7. A quadtree is a data structure in which each internal node has exactly four children.

Extension in 3D gives octree data structure where eight children exist. This data structure was named a quadtree by R. Finkel and J. L. Bentley in 1974 [START_REF] Finkel | Quad trees. A data structure for retrieval on composite keys[END_REF]. Domain discretization using square octree grids is also explained in [START_REF] Popinet | Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries[END_REF]. The tree is dened as a collection of interconnected cells, also called octants. Each renement generates disjointed subtrees, whose nodes are called children. Nodes without children are named leaves.

An octree data structure is said to be linear when only the leaves of the tree are stored in memory, which means that the usage of computational resources is optimal. While ordering the cells of a Cartesian grid is straightforward, the same cannot be said about block structured Figure 1.7: Example of quadtree discretization and corresponding tree representation. Found in [START_REF] Popinet | Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries[END_REF].

meshes, sorting is achieved using space-lling curves that cross each element of the grid once.

There are many dierent examples of space-lling curves, one of these is the Z-order, introduced by Morton [START_REF] Morton | A computed oriented geodetic data base and a new technique in le sequencing[END_REF], who also proposed a procedure, commonly known as the Morton code, to identify cells neighbours, based on simple manipulations of binary numbers. This algorithm avoids the evaluation and the storage of complex connectivity matrices, contributing to the low memory usage of octrees. The Z-order is also used to compute the parallel partition of the domain thanks to the linearity of the Z-curve. The communication between two dierent subdomains is made thanks to one layer of ghost cells where the information is shared between two processors.

Quadtree and octree allow a exibility in mesh renement and coarsening with the use of dierent levels of renement in a unique grid and also have the advantages of regular meshes as seen previously. Indeed, Cartesian grids imply a high computational eort since the whole domain has the same renement level. On the other hand other discretizations as octrees allow to rene on areas of high interest (on one body's wake, close to obstacles ...) and coarsen elsewhere. Meshing a computational domain in a complex ow context and dealing with complex geometries while making the grid t the obstacles can rapidly become very time-consuming. Moreover, changing a small part of the geometry implies doing the whole mesh again. However body-tted grids allow to get a high accuracy around obstacles which is a signicant advantage. In the context of the present work, one aim is to develop a tool that can be used by an engineer that would not be necessary a specialist in CFD. Also, for the application to the simulation of wind turbines with blades that rotate along time the use of an adaptive mesh that is not too long to compute is well adapted. Therefore, the choice to use octree grids has been done in the present work.

Existing Numerical Models in Computational Fluid Dynamics

We consider the incompressible Navier-Stokes equations in a general and bounded domain Ω that can be seen in gure 1.8. Let Ω f be the uid domain and Ω s be the solid domain, Ω = Ω f ∪ Ω s , Γ s = ∂Ω s is the boundary of the solid domain and Γ f = ∂Ω f is the boundary of the uid domain.

                   ∂u ∂t + (u • ∇)u = - 1 ρ ∇p + ν∆u in Ω f , ∇.u = 0 in Ω f , u(x, 0) = u 0 (x) in Ω f , u = u b (x, t) on Γ s , u = u f (x, t) on Γ f , (1.32
)

Ω f Ω s Γ f Γ s χ B = 1 χ B = 0 Figure 1.8: Sketch of the computational domain
where p is the pressure, u the velocity eld, ρ is the density of uid and ν the kinematic viscosity.

Spectral methods can be used in Computational Fluid Dynamics to solve Navier-Stokes equations. These methods are highly accurate in the context of incompressible ows with a simple geometry. But, when simulating with more complex geometries, a full spectral approach is not adequate. Laizet and Lamballais [START_REF] Laizet | High-order compact schemes for incompressible ows: A simple and ecient method with quasi-spectral accuracy[END_REF], developed an intermediate tool between fully spectral Navier-Stokes solvers and more versatile codes based on standard numerical schemes. They indeed used a spectral space with a Fourier representation to solve the Poisson equation, which allows to have a very accurate incompressibility condition. The other terms (convective and viscous) are discretized with high-order compact schemes. All this work is performed in the context of Cartesian grids with an immersed boundary method (this method will be discussed in the next chapter). The use of Cartesian meshes implies high computational cost when dealing with large computational domains as for the application of the ow around wind blades.

A pseudo-spectral method has also been developed by Schneider [START_REF] Schneider | Numerical simulation of the transient ow behaviour in chemical reactors using a penalisation method[END_REF]. The penalized Navier-Stokes equations in vorticity-velocity formulation are discretized in space and time. For more information about penalized equations, the reader is referred to the next chapter. The pseudospectral method is used for the space discretization, whereas for the time discretization, a variable time-stepping semi-implicit scheme is employed for the diusion term and a second-order Adams-Bashforth scheme is developed for the convective and the penalization term. The pseudospectral method consists in transforming the vorticity and other variables to Fourier space for the computation of spatial derivatives. On the other hand, the terms containing products, as the convective and penalization terms are calculated in physical space. This work has been performed in a two-dimensional case and a three-dimensional extension is required to t our present application to wind blade modelling.

The dierent discretization methods existing to solve the Navier-Stokes equations is discussed below. Schemes used for ow around wind turbines is rst investigated. Then, existing discretizations based on quadtree/octree data structure is explained. The schemes used for the advection term in these contexts are then described.

Numerical Schemes for Wind Turbine Modelling

The three-dimensional ow solver Ellipsys3D was developed in collaboration between the Technical University of Denmark [START_REF] Michelsen | Basis3da platform for development of multiblock PDE solvers[END_REF] and the former National Laboratory for Sustainable Energy (Risø) [START_REF] Sorensen | General purpose ow solver applied to ow over hills[END_REF]. This code has been developed for the purpose of modelling the wind turbine's wakes. [START_REF] Chaderjian | Advances in rotor performance and turbulent wake simulation using DES and Adaptive Mesh Renement[END_REF] the solver and its application for the modelling of helicopter's rotor which is quite close to the wind turbine application. Navier-Stokes equations are solved with central nite dierence discretizations. This code employs an overset grid approach with totally structured meshes. An overset grid system means that a near-body mesh exist in the vicinity of the rotor, body-tted and with curvilinear coordinates allowing a thin grid on the blades to resolve the boundary layer. This near-body mesh is coupled with an o-body grid system to extend the computational domain in the far eld. In the article the o-body mesh is itself composed of several Cartesian grids as "brick grids". An Adaptive Mesh Renement process (AMR) is implemented in Overow, that allows to rene the o-body grid system in areas of interest, thus, the rotor wakes can be highly resolved. An other feature of Overow is the coupling with a structural nite element model.

Commercial multi-purpose solvers as ANSYS CFX, Fluent or Star CCM + are also used in the context of ow around wind turbines. These solvers are not well documented which can be problematic in the understanding of obtained results. A multi-grid structured mesh was used in [START_REF] Laursen | 3D CFD quantication of the performance of a multi-megawatt wind turbine[END_REF] for the simulation of ow around a wind blade and nacelle, except at the rotational centre of the domain, where an unstructured grid is employed. The results obtained showed good results but some knowledge in Computational Fluid Dynamics are needed to be able to analyse the results obtained and to use adequate numerical models depending on the ow characteristics.

For all solvers described above, the meshing of the domain is not trivial and can be very timeconsuming. A particular attention has then be accorded to the discretizations of Navier-Stokes equations using an octree data structure.

Numerical Schemes with Octree Data Structure

An immersed boundary method was introduced by Peskin [START_REF] Charles S Peskin | Flow patterns around heart valves: A numerical method[END_REF] to simulate uid-structure interactions. This method with fractional-step scheme to solve incompressible Navier-Stokes equations was proposed later by Bergmann and Iollo in [START_REF] Bergmann | Modeling and simulation of sh-like swimming[END_REF] in the context of Cartesian grids. Fractionalstep methods are based on the Chorin-Temam schemes [START_REF] Chorin | Numerical solution of the NavierStokes equations[END_REF][START_REF] Témam | Sur l'approximation de la solution des équations de NavierStokes par la méthode des pas fractionnaires[END_REF]. The immersed boundary allows to avoid the use of bodytted grids, thus, this method results in a signicantly easier mesh process. Finite dierence methods are equivalent to nite volume when dealing with Cartesian grids and, in this article, space discretization is performed using a centred second order nite-dierence approximation for the diusive term. This discretization is straight-forward to implement on a Cartesian grid.

Popinet [START_REF] Popinet | Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries[END_REF] introduced a Navier-Stokes solver using an octree data structure. He proposed a second-order accurate in space and time solver on octrees. All variables are stored at the cellcentre. A Poisson equation for pressure is solved in order to maintain the second-order accuracy of the solution but the gradients have a rst-order accuracy. A large stencil is used when a level jump occurs. The variables are indeed rebuilt as ghost values using neighbours of neighbours of the cell. A constraint named grading condition is also imposed on the grid with the fact that levels of diagonally neighbouring cells can not dier by more than one.

Min and Gibou [START_REF] Min | A second order accurate projection method for the incompressible NavierStokes equations on non-graded adaptive grids[END_REF] developed a solver for Navier-Stokes equations on quadtree and octree data structures. They used the supra-convergent Poisson solver described in [START_REF] Min | A supra-convergent nite dierence scheme for the variable coecient poisson equation on non-graded grids[END_REF] with the same authors. Variables are stored at the cell's node which imply that schemes are straightforward to implement. This node-based discretization has a quite wide stencil with the use of neighbours of neighbours. It is demonstrated that second order accuracy in the L ∞ norms is obtained for the solution of Poisson equation and its gradients. The diusion term is treated with a stable backward nite dierence scheme.

In [START_REF] Olshanskii | An octree-based solver for the incompressible NavierStokes equations with enhanced stability and low dissipation[END_REF], Olshanskii et al. developed a nite dierence solver on staggered grid. The pressure variable is thus dened cell-centered and velocity component are at cell face. The advantages of using staggered location of unknowns are the enforcement of the incompressibility condition and the well-known pressure stability of such schemes. But, that implies more diculty in building higher order methods on octree grids. When a level jump occurs, the nite dierence pressure gradient uses a stencil with only the four smaller neighbours and the larger cell sharing the face.

Losasso et al. [START_REF] Losasso | Simulating water and smoke with an octree data structure[END_REF] developed an octree data structure method where pressure variable is stored at the cell centre, while the velocity components are dened on the cell faces. At the level jumps, a ghost value of pressure is computed in order to get pressure gradient at face center. A standard central dierencing is used even if it does not dene the pressure gradient exactly at cell face. It is still shown that this approximation is convergent with rst order accuracy. It has however a strong mesh constraint, since this method only allows one face with a level jump.

Another approach with a nite volume scheme was described in [START_REF] Batty | A cell-centred nite volume method for the poisson problem on non-graded quadtrees with second order accurate gradients[END_REF] by C. Batty. Batty's stencil is based on diagonals instead of axis directions as previously seen and is wide with the use of neighbours of neighbouring cells. This method allows to compute directly the ghost value at level jumps instead of doing an interpolation rst. Quadratic interpolations are constructed, which allow to get a second-order accuracy of the solution and its gradients. This method is however only described in the two-dimensional case and can not be directly applied in 3D.

Numerical Schemes for the Advection Term

The left-hand side of eq. (1.32) corresponds to the advection term of Navier-Stokes equations. Dierent schemes for its treatment exist in the literature. The semi-Lagrangian method is a widely used model. It consists of being in an Eulerian framework but with discrete equations coming from the Lagrangian perspective. It involves backward time integration of a characteristic equation to nd the departure point of a uid particle arriving at an Eulerian grid point. Interpolation formulas are then used to recover the value of solution at such points. The advantage of semi-Lagrangian compared to Lagrangian framework is the fact that a well-dened grid is conserved. Indeed, a Lagrangian method, that follows the ow along time, could result in large distortions of the mesh. Then, semi-Lagrangian methods allow to avoid computing the convective term. A feature of semi-Lagrangian schemes is the unconditional stability. As we have a time step restriction, when using explicit schemes, imposed by the CFL condition to be proportional to the smallest grid cell, the stability of semi-Lagrangian scheme removes that restriction. Several examples of semi-Lagrangian methods implemented are reported in the literature [START_REF] Xiu | A semi-Lagrangian high-order method for NavierStokes equations[END_REF][START_REF] Oliveiro | A comparison of integration and interpolation Eulerian-Lagrangian methods[END_REF][START_REF] Min | A second order accurate projection method for the incompressible NavierStokes equations on non-graded adaptive grids[END_REF][START_REF] Guittet | A stable projection method for the incompressible NavierStokes equations on arbitrary geometries and adaptive quad/octrees[END_REF][START_REF] Lentine | An unconditionally stable fully conservative semi-Lagrangian method[END_REF].

Xiu and Karniadakis [START_REF] Xiu | A semi-Lagrangian high-order method for NavierStokes equations[END_REF] have developed a second order accurate semi-Lagrangian method. Indeed, according to Oliveiro and Baptista [START_REF] Oliveiro | A comparison of integration and interpolation Eulerian-Lagrangian methods[END_REF], the simplest semi-Lagrangian scheme with linear interpolation is equivalent to the classical rst-order upwinding scheme, which is excessively dissipative. Their second order accurate scheme shows an accuracy of the backward integration as well as the accuracy of the interpolation method. An implicit midpoint rule is used for backward integration. That means that the middle point in space between the arrival and the departure point is considered. And also the mid time step t n+1/2 is employed for the characteristic curve. The velocity at the mid time step is dened as a linear combination of the velocities at the two previous time steps t n and t n-1 . The interpolation procedure to compute the departure point should then be at a second-order accuracy so that the whole semi-Lagrangian method described is itself second-order accurate.

Min and Gibou [START_REF] Min | A second order accurate projection method for the incompressible NavierStokes equations on non-graded adaptive grids[END_REF] used the aforementioned method in the context of computations of Navier-Stokes equations on an octree data structure. A similar semi-Lagrangian second-order backward scheme is developed in [START_REF] Guittet | A stable projection method for the incompressible NavierStokes equations on arbitrary geometries and adaptive quad/octrees[END_REF] for the computation of advection term on octrees. But with an adaptive time step, which requires some slight adaptation of previously described method.

Lentine et al. in [START_REF] Lentine | An unconditionally stable fully conservative semi-Lagrangian method[END_REF] noticed that semi-Lagrangian schemes are based on characteristic tracing and interpolation, thus a fully conservative implementation is not necessarily guaranteed. They proposed a novel technique that applies a conservative limiter to the typical semi-Lagrangian interpolation step in order to guarantee that the amount of the conservative quantity does not increase during this advection.

The other common method consists in an Eulerian discretization of the advection term. Popinet in [START_REF] Popinet | Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries[END_REF], used a rst-order upwind scheme with nite dierences. Olshanskii et al. [START_REF] Olshanskii | An octree-based solver for the incompressible NavierStokes equations with enhanced stability and low dissipation[END_REF], described a higher order method with the use of a third-order upwind discretization stencil with a nite dierence scheme. The results obtained show that the method used is both stable and lowdissipative.

The solver Overow, based on curvilinear or Cartesian grid systems employs discretizations of the convective term with fth-order accurate central dierences [START_REF] Chaderjian | Advances in rotor performance and turbulent wake simulation using DES and Adaptive Mesh Renement[END_REF]. The Ellipsys3D code also discretized the convective term thanks to hybrid scheme combining a third-order QUICK scheme and fourth-order CDS scheme [START_REF] Sorensen | Simulation of wind turbine wakes using the actuator line technique[END_REF]. This higher-order accuracy on both solvers can be reached without too many eorts thanks to the block structured grid approach.

Chapter 2

Numerical Modelling of Flow around a Wind Turbine on Octree Grids

As described in the previous section, quadtree (in 2D) and octree (in 3D) allow a exibility in mesh renement and coarsening with the use of dierent levels of renement on a unique grid and also have the advantages of regular meshes with a better convergence and higher resolution compared to unstructured grids. This octree-based discretization allows to rene on areas of high interest (on the body's wake, close to obstacles ...) and coarsen elsewhere, which is of particular interest for high-delity modelling of the ow around a wind turbine.

In this work, discretization and parallelization are based on a octree library named bitpit.1 More precisely, bitpit is an open source modular library for scientic computing developed by Optimad Engineering Srl. The goal of bitpit is to ease the burden of writing scientic programs providing the common building blocks needed by every scientic application. In bitpit, only the PABLO (PArallel Balanced Linear Octree) module was used. This module allows to deal with parallel linear quadtrees/octrees, thanks to the Z-order space-lling curve and dynamical adaptive mesh renement.

This module permitted an easy handling of parallel variables. Indeed, each cell surrounded by cells known by other processors (named ghost cells), has access to the so-dened ghost data (variables stored in the ghost cells). PABLO's routines allow immediate inclusion within a code without needing important adaptation. Another advantage of PABLO module is the fact that it manages linear octrees which allows a relevant gain of memory compared to other existing open-source libraries. It was also chosen in the present work that a cell has only access to the rst cells close to it (let's say the rst "layer" of surrounded cells). This implies in the context of parallelization that less communications and thus less computational time. The stencil used for the discretization is not wide and can only use the rst neighbours surrounding each cell. This chapter describes the implemented Navier-Stokes solver in the context of octree grids. The outline is, rst, the explanation of the immersed boundary method for the computation of the equations. Then, the two-dimensional nite volume solver is described, followed by the extension to the three-dimensional case. Finally, the methods developed to deal with high-Reynolds number ows are explained.

Immersed Boundary Method

Immersed boundary method is included in the larger scope of the Fictious domain method. The basic idea of ctious domain method is to substitute a given problem posed on a domain D, with a new problem posed on a simple shaped domain Ω containing D. As previously introduced, the immersed boundary method allows to have a computational grid that is not body-tted, which is simpler to implement. Indeed, as the interface between the uid and the body is included in the cell meshes, simple grids like Cartesian or octree grids can be used, which implies a low amount of memory and an easy way to parallelize, compared to unstructured body-tted meshes. Two approaches exist in the immersed boundary method: continuous forcing and discrete forcing. In the former, a force term is added to the continuous Navier-Stokes equations before discretization, whereas in the latter, the forcing is applied (explicitly or implicitly) to the discretized equations. The main advantage of the discrete forcing, such as in the ghost-cell method proposed by Mittal et al. [START_REF] Mittal | A versatile sharp interface immersed boundary method for incompressible ows with complex boundaries[END_REF], is the possibility to achieve a sharp representation of the uid-solid interface, exactly as if the mesh were body tted. On the other hand, its drawback is the problem of the so-called "fresh cells", a situation that is encountered dealing with moving interfaces: some solid cells might emerge into the uid between one time step and another as a result of the boundary motion. Continuous forcing methods are not sharp, indeed the immersed boundary is diused, leading to a loss of accuracy in the proximity of the body, but they allow to bypass the special treatment of fresh cells, as remarked by Bergmann et al. [START_REF] Bergmann | An accurate Cartesian method for incompressible ows with moving boundaries[END_REF]. For the present work, a continuous forcing model was employed, which corresponds to the volume penalization method described in [3].

Penalized Equations

Using the notation of gure 1.8, the previously described Navier-Stokes equations (1.32) can be written using the penalization method, where the solid body is considered as a porous medium with very low permeability κ:

             ∂u ∂t + (u • ∇)u = - 1 ρ ∇p + ν∆u + χ B κ (u B -u) in Ω, ∇.u = 0 in Ω, u(x, 0) = u 0 (x) in Ω, u = u f (x, t) on Γ f , (2.1)
where χ B is the characteristic function:

χ B = 0 in the uid domain 1 in the solid domain (2.2)
It has been demonstrated by Angot et al. [3], that the system (2.1) converges to the system (1.32) when κ → 0. More precisely, Nguyen van yen et al. [START_REF] Nguyen Van Yen | Approximation of the Laplace and Stokes operators with Dirichlet boundary conditions through volume penalization: a spectral viewpoint[END_REF] proved that the convergence of the model is in √ κ.

In the present work, the solid body geometry is represented by Lagrangian markers. Basically, the body is divided in many sections and each section is formed by many points. Thanks to these markers, a signed distance function is computed for each cell of the grid. Indeed, for each point of our regular octree grid, the closest marker is searched in order to determine the distance from the surface. The Lagrangian markers have been organized, so that the sign of the distance function can be computed with the sign of the scalar product between the outward 1. Immersed Boundary Method surface normal and the vector that links the surface to the point in the domain. In order to simplify the notation, the signed distance function is called in present work the level-set. The relationship between the level-set function ϕ and the characteristic function is the following:

χ B = 1 -H(ϕ), (2.3)
with H(ϕ), the Heaviside function.

In the thesis, a rst-order penalization is used in the solver developed, as a rst approach. By the use of the level-set function gradients, the outward normal to the body is also computed at the uid-solid interface:

n = ∇ϕ ∇ϕ ϕ=0 . (2.4)
The penalized Navier-Stokes equations are then solved numerically using a fractional-step method as will be explained below.

Temporal Discretization

A semi-Lagrangian scheme has been considered in the present work which means that a Lagrangian perspective is used in the Eulerian framework. The left hand side term of the momentum equation in eq. ( 2.1) can be written:

∂u ∂t + (u • ∇)u = Du Dt . (2.5)
As we have a decoupling between the pressure and the velocity elds, nearly all numerical methods to solve the Navier-Stokes equations (2.1) use a fractional-step approach. This approach, referred to as the projection method, has been pioneered by Chorin [START_REF] Chorin | Numerical solution of the NavierStokes equations[END_REF] and Temam [START_REF] Témam | Sur l'approximation de la solution des équations de NavierStokes par la méthode des pas fractionnaires[END_REF]. For a complete review of the dierent projection methods, the reader can refer to the article of Brown, Cortez and Minion [START_REF] Brown | Accurate projection methods for the incompressible NavierStokes equations[END_REF]. In the present work, three fractional steps are solved with the computation of two intermediate velocity elds u * and u * * . The subscript a means the arrival point, that ts with the mesh, whereas the subscript d corresponds to the departure point. The momentum equation can indeed be written:

u n+1 a -u n d -u * a + u * a + u * * a -u * * a ∆t = - 1 ρ ∇p n+1 + χ B κ (u n+1 B -u n+1 a ) + ν∆u * a . (2.6)
An intermediate velocity u * is then computed as a prediction value of the velocity eld u. This rst step, named the prediction step, consists in solving the following linear system:

u * a -u n d ∆t = ν∆u * a .
(2.7)

The diusive term of the Navier-Stokes equations is implicitly computed. u n d represents the velocity eld interpolated at the root of the characteristic curve. This interpolation will be detailed in the next section. For easy handling, no initial guess of the pressure eld is included for this prediction step. This algorithm is said to be non-incremental and is more precise, since there is no cumulative error on the pressure eld. However, it is less stable than the incremental method, due to bigger correction imposed by the pressure. As far as we tested our code, no issue related to our non-incremental scheme occurred. Coming back to the prediction step, the computed predicted velocity eld u * is not supposed to satisfy the divergence constraint. The correction step consists in solving a second fractional step:

u * * a -u * a ∆t = - 1 ρ ∇p n+1 . (2.8)
Beforehand, in order to enforce the divergence constraint, the divergence operator is applied to eq. (2.8) and we obtain the following elliptic equation referred to as the projection step:

∇ • u * a = ∆t ρ ∆p n+1 .
(2.9)

Computing this equation allows to determine the pressure eld. Then, we come back to eq. (2.8) (correction step) and the computation of the velocity eld can be done. In order to apply the penalization term, a third fractional step is performed:

u n+1 a = u * * a + ∆t χ B κ (u n+1 B -u n+1 a
).

(2.10)

Boundary Conditions

Boundary conditions for the pressure and the velocity eld need to be imposed to correctly solve the problem. Considering the pressure, Neumann boundary conditions on all edges of the domain borders t our application to study the ow past an obstacle. However, that leads to a singular linear system for the Poisson equation (2.9). It was thus chosen to remove the null space (i.e. the solution with a null average) to solve the system with Neumann boundary conditions. Indeed, the equations can be solved up to a constant. Another possibility would have been to

x the pressure at one edge of the border.

The boundary conditions for the velocity eld will now be described. The inlet boundary condition will consist in imposing an incoming ow. So, Dirichlet boundary conditions are used at inlet corresponding to wind speed for our application to wind turbine modelling. On the contrary, at the outlet, Neumann boundary conditions are imposed to allow the ow to get out of the domain. The lateral borders should not let other ow getting in on the normal direction, therefore, a zero velocity Dirichlet condition is imposed for the normal component of the velocity eld. For the tangential components (one component in two-dimensional case), it was chosen to impose a zero gradient condition. The initial condition of the velocity eld is chosen to be the same as the inlet boundary condition. The Neumann boundary conditions imposed at the outlet constitute a rst approach for simplicity, but numerical results can be improved by imposing non-reecting boundary conditions instead. Especially, the size of the computational domain could be reduced.

Space and time discretizations of the Navier-Stokes solver are now described in detail for the two-dimensional case before extending the method to the three-dimensional case.

2 Two-Dimensional Finite Volume Solver A Navier-Stokes solver with immersed boundaries in the framework of octree grids has been developed during this thesis. As explained in the previous section, a fractional-step method is used and the code solves equations (2.7), (2.9) and (2.8). The prediction and projection steps consist in solving linear systems. The open-source PETSc library is used to solve the linear systems [4]. Krylov subspace methods are employed, since these methods are known to be very successful in numerical linear algebra. Since we are using a discretization on octree grids with the Z-order space-lling curves, the variables in two neighbouring cells can be located far from each other in the matrices of the linear systems. Thus, two dierent iterative methods, that are part of Krylov subspace theory, have been tested for our Navier-Stokes solver: biconjugate gradient stabilized and generalized minimal residual methods. Both gave us the same results when compared with an analytical solution. This will be explained in more details in the next chapter dedicated to the validation of the solver. The nite volume method has been widely used in the present work for the octree-based discretization of the Navier-Stokes problem since it allows an easy formulation based on the evaluation of uxes at adjacent cell surfaces. All the pressure and velocity variables are here only stored at cell-centres, since the PABLO module, used to handle parallel discretization, is more adapted to this conguration.

The prediction step (2.7) consists in solving two linear systems (for the two-dimensional case). A semi-Lagrangian scheme allows to compute the advection term. The left-hand side is constituted by a matrix constructed by the Laplace operator for the discretization of the implicit diusive term multiplied by the vector of unknowns constituted with the predicted velocity eld u * . The projection step (2.9) is an elliptic equation to obtain the pressure eld. A linear system should thus be solved with a right-hand side corresponding to the discretization of the divergence of u * . The other side of the equation is similar as previously seen since a matrix is constructed that corresponds to the implicit discretization of the Laplace operator for pressure. Finally, the correction step (2.8) is quite simple since u * and p are known. The gradient of the pressure is computed by the use of the nite volume method as will be explained below and the velocity eld u n+1 can be obtained.

Numerical Scheme for the Divergence Operator

For the correction step, the computation of the pressure gradient is needed at the cell-centre in order to get the velocity eld u n+1 . In a similar way, the divergence of u * should also be computed in the projection step. The explanation of the method used will be described for the case of the divergence of some vector V (totally similar to the divergence of u * ) for simplicity.

Exactly the same method is employed for the gradient computation. As explained above, a nite volume method is used, which means that we integrate ∇ • V over each cell volume Ω c . The divergence theorem and its corollaries are then employed to convert the volume integral into surface integrals:

1 Ω c Ωc (∇ • V)dΩ c = 1 Ω c ∂Ωc (V • n)dS, (2.11)
where n represents the outward normal of the cell. The right-hand side of eq. (2.11) is then evaluated as the sum of all uxes at the edges of each cell. Thus, we can write for each cell i of the domain:

(∇ • V) i = 1 ∆Ω i j V j .n j ∆S j , (2.12)
with j representing the edge index. For the two-dimensional case ∆S j corresponds to the edge length and ∆Ω i to the area of i th cell. V j needs to be evaluated for each cell at the edge centre.

In order to optimize the operations, the developed solver computes only once the ux at each edge and uses this ux for the two surrounded cells sharing the edge.

When no level dierences exist between the two cells sharing the edge (so when dealing with a Cartesian case), V j .n j is basically computed as the average of the two neighbouring cells where V is known at cell-centre. In case of a level jump, for example between the cells c 4 and c 1 as can be seen in gure 2.1, V j .n j is evaluated with the vertex belonging to the edge:

V j .n j = V v 1 + V v 2 2 at level jump. (2.13)
But V v 1 and V v 2 are unknowns at cell vertices, thus, their values need to be interpolated. All the neighbours surrounding each vertex are searched and linear or bilinear interpolation is then

+ + + + + o x x n c 1 c 2 c 3 c 4 c 5 v 1 v 2 Figure 2
.1: Example of mesh conguration in 2D with stencils computed with the available neighbours. Indeed, only two dierent congurations exist: either the vertex v 1 has 4 neighbours which implies that vertex v 2 has 3 neighbours, or the opposite. A special treatment is done when the level jump has a vertex located at the border of the domain. In this case, the interpolation is computed at the edge centre with a stencil composed of the neighbours sharing the edge. In the present case with a level jump, we can write:

V v 1 = α v 1 + β v 1 .x v 1 + γ v 1 .y v 1 + ζ v 1 .x v 1 .y v 1 , V v 2 = α v 2 + β v 2 .x v 2 + γ v 2 .y v 2 , (2.14)
with x v i and y v i that represent the coordinates of the vertex v i . Small linear systems are solved to obtain the unknown coecients

α v i , β v i , γ v i and ζ v i .        1 x c 1 y c 1 x c 1 .y c 1 1 x c 2 y c 2 x c 2 .y c 2 1 x c 3 y c 3 x c 3 .y c 3 1 x c 4 y c 4 x c 4 .y c 4               α v 1 β v 1 γ v 1 ζ v 1        =        V c 1 V c 2 V c 3 V c 4        , (2.15)
    1 x c 1 y c 1 1 x c 4 y c 4 1 x c 5 y c 5         α v 2 β v 2 γ v 2     =     V c 1 V c 4 V c 5     .
(2.16)

In order to proof the consistency and to evaluate the order of this method, an analysis has been carried out. The function f (x, y) = sin(x 2 + y 2 ) has been chosen in the two-dimensional square domain Ω = [-0.5, 0.5] × [-0.5, 0.5]. Exact Dirichlet boundary conditions have been imposed at the cell centres along the border of the domain. The computation of ∇.f has so been performed by the use of the previously described solver and has been compared to the exact solution of the divergence of f . This analysis has been made for a Cartesian grid as a rst validation and the results can be seen in table 2.1. As expected, an order of convergence of 2 is obtained for both the L 2 and the L ∞ norms. The same study has then been performed with a quadtree grid. This grid has been chosen with a high number of level jumps in order to prove the consistency The results can be seen in table 2.2. Thus, the L ∞ norm is consistent with an order of 1. L 2 norm has a 1.5 order accuracy. We even get the second-order accuracy when looking at the L 1 norm of the error. This divergence operator scheme is so consistent even with highly irregular quadtree grid conguration.

+ + + + + o x x n c 1 c 2 c 3 c 4 c 5 v 1 v 2 τ c τ v

Numerical Scheme for the Laplace Operator

The linear systems that should be solved for both the prediction step (2.7) and the Poisson equation (2.9) imply the discretization of a Laplace operator. To explain the method developed in the present work, the Laplacian of a variable Φ will be described for easier writing. The reader is reminded that exactly the same method is used for the discretization of the pressure Laplacian in (2.9) and of the predicted velocity eld Laplacian in (2.7). A nite volume method is employed and we integrate ∇ 2 Φ on each cell volume. The divergence theorem is then applied:

1 Ω c Ωc ∇ • (∇Φ)dΩ c = 1 Ω c ∂Ωc (∇Φ • n)dS.
(2.17)

The right-hand side of eq. (2.17) is then evaluated as the sum of uxes for each edge of cells, which can be written for each cell i of the domain:

∇ 2 Φ i = 1 ∆Ω i j (∇Φ) j .n j ∆S j . (2.18)
The problem comes down to the evaluation of the gradient of Φ at each cell edge. Discrete Duality Finite Volume (DDFV) schemes have been investigated. The aim of such methods is therefore to provide a discrete reconstruction of gradients. The idea of these methods is to integrate the equations on both a primal and a dual mesh. Coudière et al. [START_REF] Coudière | Convergence rate of a nite volume scheme for a two dimensional convection-diusion problem[END_REF] developed the reconstruction of discrete gradients on a diamond mesh and this method has also been implemented on the present work since the mesh conguration with quadtrees could t with diamond dual mesh. An example of grid can be seen in gure 2.3. When having no level dierence between the two cells sharing the edge, for example, when computing the gradient on the edge shared by c 3 and c 2 , the gradient can simply be computed:

∇Φ j .n j = Φ c 3 -Φ c 2 ∆l
• n j for a Cartesian conguration (2.19) with ∆l corresponding to the length of the edge.

When dealing with a quadtree conguration, a diamond cell is used in order to reconstruct the discrete gradient of this dual mesh. The following system should so be solved:

     ∇Φ. τ c = Φ c 1 -Φ c 4 ∆c , ∇Φ. τ v = Φ v 1 -Φ v 2 ∆l , (2.20)
with ∆c being the length between the centres of cells c 1 and c 4 . It can then be written that:

∇Φ = ∇Φ 1 . e 1 + ∇Φ 2 . e 2 = ∇Φ 1 . n + ( ∇Φ. τ v ). e 2 , (2.21)
with e i that represents the two unit vectors in x or y direction. When dealing with quadtree grids, τ v is either oriented in x or y direction and n in the other direction. In the present case we chose to write e 1 = n. That implies also that:

∇Φ. n = ∇Φ 1 .
Thanks to the use of (2.20) and (2.21), we nally obtain:

(∇Φ) j .n j = 1 τ c .n j Φ c 1 -Φ c 4 ∆c - Φ v 1 -Φ v 2 ∆l . τ v . τ c at level jump. (2.22)
In the present work, it was chosen to store only variables at the cell-centre. Moreover, we are dealing with an implicit variable. So, Φ v 1 and Φ v 2 should be evaluated from cell-centred vari- ables. In order to use the same stencil as the previously described divergence operator, linear or bilinear interpolations were also employed for the present case. The same linear systems as eqs (2.15) and (2.16) can thus be written for the variable Φ but they can't be directly solved since we are dealing with implicit variables. The unknown coecients are expressed from Φ variable and the matrix corresponding to the global linear system of the prediction or projection steps can then be constructed.

This scheme for discretization of the Laplace operator has then been analysed to evaluate its consistency and accuracy. The same analytical function, same square domain and the same grid congurations as in the previous section have been used and the equation ∇ 2 f = s is solved with s(x, y) being the analytical Laplacian of f (x, y) = sin(x 2 + y 2 ). Exact Dirichlet boundary conditions have been imposed at the cell centres along the border of the domain. A rst analysis has been carried out with a Cartesian grid to proof the consistency of the method. The results can be seen on table 2.3. We reach an order of convergence of 2 with both the L 2 and L ∞ norms, which is the expected result.

The same study has also been made with the quadtree grid conguration previously used (gure 2.2). The results obtained can be observed on table 2.4. This analysis shows a second-order accuracy with both the L 2 and L ∞ norms for the Laplace operator discretization even on a highly irregular grid conguration. However, the level of errors are bigger compared to the Cartesian conguration. That is explained with the fact that no interpolation is needed for the Cartesian grid compared to quadtree. The results obtained with the L 1 norm are not presented here but similar results as for the L 2 norm are noticed. Coudière et al. [START_REF] Coudière | Convergence rate of a nite volume scheme for a two dimensional convection-diusion problem[END_REF] explained however that on general meshes, an order of convergence of 1 is reached for the Diamond's method. Thus, we might have a superconvergence phenomenon with the Laplace operator discretization. An example of the error distribution on a rened grid is presented in gure 2.4. The exact solution is symmetric in the computational domain but the error distribution is not. This can be explained by the fact that the grid conguration is not symmetric in that domain. The Laplace operator scheme is thus consistent and shows a second-order accuracy even with highly irregular grid conguration.

+ + + + + + + + x c 1 c 2 c 3 c 4 c 5 c 6 c 7 c 8 u n (x c1 -∆t.u n (x c1 ))
x d

Semi-Lagrangian Scheme for the Advection Term

The Lagrangian description of a ow eld requires to follow the trajectories of the ow particles whereas the Eulerian description xes a particular point in space and considers the time evolution of variables at this point. In this work, a semi-Lagrangian description of the ow is adopted. It avoids the direct computation of the convective terms which is not straightforward with an Eulerian scheme in the presence of non-conforming meshes as quadtree grids. The semi-Lagrangian scheme consists in computing the predicted velocity eld by taking into account the diusive term and the velocity computed at the previous time step by following the trajectory of the eld. So, a semi-Lagrangian scheme uses an Eulerian framework but the discrete equations come from the Lagrangian perspective. Thus, in the prediction step, the characteristic curve is computed for each cell i as:

(x n d ) i = x i -∆tu n i , (2.23)
with x i being the position of i th cell centre, u n i corresponds to the velocity eld at cell centre and (x n d ) i represents the position at the root of the characteristic. This quantity needs to be interpolated since there is no reason that this position coincides with the grid. As this characteristic curve uses only a single time-step backward integration, the time discretization of the Navier-Stokes solver will be only rst-order accurate in time.

As previously explained, by the use of the bitpit library, we have a constraint in parallel computing, which is the fact that each cell knows only its rst neighbouring cells. Indeed, the implicit computation of the diusive term in the Navier-Stokes equations (2.1) prevents us to be constraint by the stability condition and allows us to use large time steps. Moreover the semi-Lagrangian scheme is unconditionally stable. However, the parallel computing imposes that the root of the characteristic stands at least in a ghost cell in order to interpolate the velocity.

The method developed for the implementation of the semi-Lagrangian scheme in the solver will now be described. For more simplicity the description will be done with the grid conguration shown in gure 2.5. The idea to interpolate the velocity eld at the characteristic root is to have a loop on each cell, then obtaining the index of the cell where to interpolate and nally having a plane reconstruction of the velocity eld in that cell. In the present conguration, our loop is currently on cell c 1 and the characteristic curve root is located on cell c 4 . So, a plane equation is used to rebuild the velocity eld in the cell:

u n • e i = α + βx + γy, (2.24)
with u.e i representing the two components of the velocity eld. It will be simplied by looking at one component. The three parameters are then found thanks to the use of the velocity at the cell-centre (which is known) and with the gradient computation at cell-centre. We indeed have:

         ∂u ∂x c 4 = β, ∂u ∂y c 4 = γ, (2.25)
and

u n c 4 = α + ∂u ∂x c 4 .x c 4 + ∂u ∂y c 4 .y c 4 . (2.26)
Finally, we obtain:

u n (x) = u n c 4 - ∂u ∂x c 4 .x c 4 - ∂u ∂y c 4 .y c 4 + ∂u ∂x c 4 .x + ∂u ∂y c 4 .y. (2.27)
The computation of the gradients at each cell-centre is done with a nite volume method. Exactly the same method as previously described for the divergence operator is used. This computation is done previously to the loop for semi-Lagrangian and the gradient values are stored. This scheme to treat the advection term in the Navier-Stokes solver has been tested with an analytical function which is transported in a quadtree grid conguration. The aim is to evaluate the interpolation error at the root of the characteristic curve. The same meshes have been used and the analytical function is g(x, y, t) = f (x, y).e t with t being the time and f (x, y) the previously described sinus function. The results of errors and order of convergence for the space discretization of the semi-Lagrangian can be observed on table 2.5. A rst-order accuracy is reached for the L ∞ norm for a highly non-conformal grid and we reach even a second-order accuracy with other norms. Thus, we validated all the operators needed to solve the whole Navier-Stokes equations for the two-dimensional case. In the next chapter, the tool has been widely tested with several two-dimensional benchmark cases existing in the literature.

3 Extension of the Methods to the Three-Dimensional Case

The methods previously described have been chosen in order that the extension to the threedimensional case can be straight-forward. In the same way as for the two-dimensional solver, the dierent operators have been developed separately and validated thanks to the comparison with an analytical function before putting all together to solve the whole Navier-Stokes equations (2.1). The main change done in the three-dimensional solver compared to the two-dimensional one is the way we deal with interpolations. Indeed, polynomial interpolations in 3D can be painstaking because of the large number of possible congurations that can be created depending of the number of neighbours. Therefore, Gaussian Radial Basis Function (RBF) interpolation has been chosen:

f (x) = i w i e r 2 i , (2.28)
where w i are the weights of the interpolation, r i = xx i and is a scaling parameter proportional to 1 ∆l i (inverse of the cell sizes). This method is second-order accurate and is independent on the grid conguration, which implies an easy implementation. One of the drawbacks is the fact that it relies only on the distance between the desired point and the interpolation points. The position is not taken into account. However, our tests showed satisfying results with this method.

Numerical Scheme for the Divergence Operator

The method used to discretize the divergence operator in 3D is the same as for the twodimensional case. That means that a nite volume method is used (see eq. (2.11)). We can write for each cell i of the domain:

(∇ • V) i = 1 ∆Ω i j V j .n j ∆S j , (2.29)
with j representing the face index. For the three-dimensional case ∆S j corresponds to the face area and ∆Ω i to the volume of the i th cell. V j needs to be evaluated at face centre for each cell. When no level dierence occurs between the two cells surrounding the face, V j .n j is basically computed as the average of the neighbouring cells where V is known at the cell-centre. In case of octree conguration, the rst idea tested was based on the two-dimensional method (see eq. (2.13)) that consists in, rstly, interpolating on the vertices composing the face (with RBF-based interpolations) and secondly, coming back to face-centre interpolation. However, when tested, this method showed inaccurate results. This was perhaps due to the fact that RBF interpolations don't take into account the position of the point used in the stencil. So, the method developed in the solver is a direct interpolation of the variable at face-centre with a single RBF interpolation. The stencil used is quite wide since all the neighbours of each vertex constituting the face are used. A criterion is however added to retrieve from the stencil the neighbours that are too far from the face-centre (> 1.5∆l with ∆l the length of the face). It was indeed noticed with the use of RBF interpolations that a stencil too large implied less accurate results. This operator has been analysed in order to proof the consistency and to evaluate the order of convergence with irregular octrees. The mesh congurations used can be seen in gure 2.6. The function f (x, y, z) = sin(x 2 + y 2 + z 2 ) has been chosen in the three-dimensional domain Ω = [-0.5, 0.5]×[-0.5, 0.5]×[-0.5, 0.5]. Exact Dirichlet boundary conditions have been imposed at the cell centres along the borders of the domain. On a Cartesian grid, the error norms show second-order accuracy. The results obtained with the octree grids can be seen in table 2.6. It Table 2.6: Norm of the error and order of the divergence operator for 3D octree grid can be observed that the L ∞ norm of the error is rst-order accurate. L 2 norm has a 1.5 order of accuracy and a second-order accuracy is obtained for the L 1 norm. This allows to prove the consistency of the divergence operator scheme with highly irregular octree grid congurations.

Numerical Scheme for the Laplace Operator

The linear systems that should be solved for both the prediction step (2.7) and the elliptic equation (2.9) imply the discretization of a Laplace operator. An extension of the method used for the two-dimensional case has been done and a nite volume scheme is used (see eqs.

(2.17) and (2.18)). The gradients need so to be evaluated at each face centre. For a Cartesian conguration, exactly the same method as in 2D case is used with the use of the two surrounding neighbours (see eq (2.19)). When dealing with an octree conguration, a diamond cell is used to reconstruct the discrete gradient of this dual mesh at the centre of each face. An example of

+ + o c 1 c 0 v 0 v 1 v 2 v 3 τ 1 τ 2 Figure 2
.7: Example of mesh conguration in 3D with stencils for Laplace operator the grid can be seen in gure 2.7. The following system should be solved:

                 ∇Φ. τ c = Φ c 0 -Φ c 1 ∆ c , ∇Φ. τ 1 = 1 2 ∇Φ. τ v 3 v 1 + ∇Φ. τ v 2 v 0 = 1 2 Φ v 3 -Φ v 1 ∆l + Φ v 2 -Φ v 0 ∆l , ∇Φ. τ 2 = 1 2 ∇Φ. τ v 3 v 2 + ∇Φ. τ v 1 v 0 = 1 2 Φ v 3 -Φ v 2 ∆l + Φ v 1 -Φ v 0 ∆l , (2.30)
with ∆c being the length between the cell centres c 0 and c 1 , and ∆l the size of the face.

Using the same methodology as for the two-dimensional case, we can nally write:

(∇Φ) j .n j = 1 τ c .n j Φ c 0 -Φ c 1 ∆c - τ 1 . τ c 2∆l (Φ v 3 -Φ v 1 + Φ v 2 -Φ v 0 ) - τ 2 . τ c 2∆l (Φ v 3 -Φ v 2 + Φ v 1 -Φ v 0 ) .
(2.31) The variables are only known at cell centres, so, the variables at nodes should be interpolated to apply this scheme. All the neighbours of each vertex are used in the stencil and a RBF interpolation is performed at each node with the method described previously. The consistency and accuracy of this scheme in 3D has also been evaluated by solving a heat equation with an analytical function: ∇ 2 f = s with s(x, y, z) being the analytical Laplacian of the previously used function f (x, y, z) = sin(x 2 + y 2 + z 2 ). Exact Dirichlet boundary conditions are imposed at the borders of the domain Ω. The results on error can be observed on table 2.7. The orders of convergence with both the L 2 and L ∞ norms are second-order accurate. This allows to validate our Laplacian operator even on highly irregular octree grids.

Semi-Lagrangian Method

For the three-dimensional solver, a semi-Lagrangian scheme has been chosen to deal with the advection term, as for the 2D code (see eq. (2.23)). The velocity needs to be interpolated at the root of the characterisitc curve since there is no reason that this point coincides with the grid. For the interpolation it was chosen to perform RBF interpolations as it seemed to be the easier Table 2.8: Norm of the error and order of the semi-Lagrangian scheme for 3D octree grid method. Indeed, the closest vertex to the root of the characteristic curve is obtained and all the cells surrounding this vertex are used in the stencil. This scheme has then be tested with an analytical function which is transported in a octree mesh conguration. The same meshes are used and the analytical function is g(x, y, z, t) = f (x, y, z).e t with t being the time and f (x, y, z) the previously described sinusoidal function. The results of errors and order of convergence of the space discretization can be seen on table 2.8. A second-order accuracy is reached for all norms. A CF L condition of 0.7 has been chosen since it allowed to obtain a good compromise between accuracy and future time computing.

The time discretization of our 3D solver is however only rst-order accurate considering the semi-Lagrangian scheme used as explained by Falcone and Ferretti [START_REF] Falcone | Convergence analysis for a class of high-order semi-lagrangian advection schemes[END_REF]. Finally, all the operators needed for solving the Navier-Stokes in 3D have been validated. The whole solver has then been tested with several benchmark cases existing in the literature as will be explained in the next chapter.

Wall Functions for High-Reynolds Number Flow

Characteristic Reynolds numbers of a wind turbine ow in operative conditions are of the order of millions [START_REF] Ceyhan | Towards 20MW wind turbine: High Reynolds number eects on rotor design[END_REF]. That means that the uid motion is characterized by the presence of turbulent boundary layers and multiscale vortical structures in the wake. Thus, performing a direct numerical simulation won't be possible. Given the highly unsteady nature of the ow that will be simulated, a promising approach is the use of Large Eddy Simulation (LES) that permits to resolve only the largest structures of the turbulence. Indeed they carry the largest part of the energy and give the biggest contribution to the transport phenomena. The smallest structures, that are associated with turbulent energy dissipation are modelled using a sub-grid model. That method will now be explained.

Large Eddy Simulations

Based on the book of S. Pope [START_REF] Stephen | Turbulent Flows[END_REF], the mathematical formulation of LES is based on a low pass lter operator (noted G), so that ltered variables are dened as:

ψ = x 1 x 2 x 3 G(r, x)ψ(x)dx 1 dx 2 dx 3 .
(2.32)

Applying (2.32) to Navier-Stokes equations (1.32) allows to obtain ltered equations. Considering a spatially uniform lter, ltering operation commutes with dierentiation. The continuity equation is linear, therefore ltering the equation is equivalent to apply the divergence operator to the ltered variables:

∂u i ∂x i = ∂u i ∂x i = 0, (2.33)
where Einstein notation is applied. The ltered momentum equation is:

du j dt + ∂u i u j ∂x i = - 1 ρ ∂p ∂x j + ν ∂ 2 u j ∂x i ∂x i . (2.34)
Due to the non-linearity of the convective term, the form of the equation is dierent from the non-ltered one. Indeed, the ltered product U i U j is not equal to the product of ltered velocities. The residual-stress tensor is dened as the dierence of the two elements:

τ ij = u i u j -u i u j , (2.35)
so that equation (2.34) can be rewritten as:

du j dt + ∂u i u j ∂x i = - 1 ρ ∂p ∂x j + ν ∂ 2 u j ∂x i ∂x i + ∂τ ij ∂x j . (2.36)
From a physical point of view, the residual-stress tensor can be interpreted as the exchange of momentum, at the ltered scale, exerted by the sub-grid turbulent structures. Let k r be the residual kinetic energy:

k r = 1 2 τ ii . (2.37)
Then, the residual-stress tensor can be decomposed into an isotropic component and an anisotropic one:

τ r ij = τ ij - 2 3 k r δ ij , (2.38)
where δ ij represents the Kronecker delta. The isotropic component is included in the modied ltered pressure term:

p * = p + 2 3 ρk r . (2.39)
The ltered momentum equation can now be rewritten into its nal form:

∂u j ∂t + ∂u i u j ∂x i = - 1 ρ ∂p * ∂x j + ν ∂ 2 u j ∂x i ∂x i + ∂τ r ij ∂x j . (2.40)
In order to solve the equation, the problem needs to be closed using what is called a subgrid model, which express the anisotropic residual-stress tensor as a function of the ltered variables. Most of the sub-grid models rely on the Boussinesq hypothesis, which states that τ r ij is proportional to the ltered rate of strain tensor:

τ r ij = ν e S ij = ν e ∂u i ∂x j + ∂u j ∂x i , (2.41)
where ν e is the eddy viscosity. The rst and the simplest sub-grid model is the Smagorinsky model [START_REF] Smagorinsky | General circulation experiments with the primitive equations[END_REF]:

ν e = (C S ∆) 2 S, (2.42)
where C S is the Smagorinsky constant, whose classical value is 0.17, ∆ is the length of the lter, with:

S = 2S ij S ij . (2.43)
This model is thus very simple to implement, but its main drawback is the poor performance in transitional ows. Indeed, the appropriate value of C S depends on the ow regime. Many eorts have been made to extend this model and one of the most notable is the dynamic model proposed by Germano et al. [START_REF] Germano | A dynamic subgrid-scale eddy viscosity model[END_REF]. This model is able to predict with accuracy regions of laminar, transitional and fully developed turbulent ow, but it is not easy to implement for complex threedimensional geometries. For the present work, it was chosen to employ the Vreman model [START_REF] Vreman | An eddy-viscosity subgrid-scale model for turbulent shear ow: Algebraic theory and applications[END_REF]:

ν e = c B β α ij α ij , (2.44)
with:

α ij = ∂u j ∂x i , (2.45)
β ij = ∆ 2 m α mi α mj , (2.46)
B β = β 11 β 22 -β 2 12 + β 11 β 33 -β 2 13 + β 22 β 33 -β 2 23 , (2.47)
and c 2.5C 2 S .

(2.48)

This model is simpler to implement than the dynamic one and it is equally capable of predicting the various characteristics of dierent ow regimes.

The ltered system of equations is thus closed by applying the Boussinesq hypothesis to the anisotropic residual-stress tensor, whose contribution to the momentum equation is:

∂τ r ij ∂x j = ∂ ∂x j ν e ∂u i ∂x j + ∂u j ∂x i . (2.49)
Using the Vreman model, the eddy viscosity is a function of ltered velocities. Hence, the turbulent stress term is decomposed into an additional viscous term and two highly non-linear terms related to the gradient of the eddy viscosity:

∂τ r ij ∂x j = ∂ν e ∂x j ∂u i ∂x j + ∂ν e ∂x j ∂u j ∂x i + ∂ 2 u j ∂x i ∂x i . (2.50)
These contributions of the gradient of the eddy viscosity are usually very small and they are therefore neglected in the global computation.

The model has been implemented in the three-dimensional solver as an additional prediction step:

u * * i = u * i + ∆t ∆Ω i ν e LU ∞ j ∂u * ∂n j ∆S j . (2.51)
The model comes down to compute discrete gradients of velocity eld at each face centre. The diamond's method, as previously explained, has been used. In the simulations performed in the present work, the classical value of the Smagorinsky constant, equal to 0.17 has been used. that the number of grid points required by a wall resolved LES is of the order of Re 13 7 , whereas it is approximately equal to Re for wall modelled LES. Moreover, the use of octree discretization in the developed solver, prevents to properly model the boundary layer near walls. Indeed, grid renement with octrees allows only to divide by 2 the size of the cells whereas the anisotropic boundary layer around the wind turbine blades would require to have very thin cells around the obstacle. Wall modelling has thus been investigated in the present work. The simplest approach for wall modelling is the use of the so-called wall functions, by which it is possible to impose a correction on the predicted velocity eld in the rst cells close to the uid-body interface. The idea behind this technique is that the turbulent boundary layer velocity prole on a at plate and with no external pressure gradient, is universal if properly normalized. Let u τ be the friction velocity (or shear velocity), dened as:

u τ = τ w ρ , (2.52)
with τ w the wall shear stress, that express:

τ w = µ ∂U ∂y w , (2.53)
where U is the tangential velocity and y is the normal direction. Using the friction velocity it is possible to dene two adimensional variables by which an universal velocity prole can be described: The outer layer is composed by part of the logarithmic layer and the wake region and the ow is there strongly dependant on the geometry and on the Reynolds number. On the other hand, the inner layer, formed by the viscous sublayer and the buer layer, allows to dene a function f so that U + = f (y + ). Traditionally, the universal prole has been described using two distinct analytical functions:

y + = y u τ ν , U + = U u τ . ( 2
u + = f (y + ) =    y + , if y + < 5 1 k logy + + B, if y + > 30, y δ < 0.3 (2.55)
where k = 0.41 is the Von Karman constant, B = 5.2 and δ is the boundary layer thickness, which is a function of the Reynolds number. The buer layer is a transition region between the two laws. It is possible to nd in the literature many dierent functions to describe the prole with a unied law, such as those developed by Spalding [START_REF] Spalding | A single formula for the Law of the Wall[END_REF]. For the present work, the Reichardt's wall law is considered as described in [2]:

u + = f w (y + ) = 1 k log(1 + ky + ) + 7.8 1 -e y + 11 - y + 11 e -0.33y + . (2.56)
As it can be noticed in gure 2.9, the advantage of using Reichardt's wall law is the possibility to capture with a reasonable level of accuracy the three dierent regions of the inner layer using one analytical expression.

The numerical implementation has been done in the three-dimensional solver by following the guideline of De Tullio [START_REF] De Tullio | Development of an immersed boundary method for the solution of the preconditioned NavierStokes equations[END_REF]. In a rst loop, all the interface cells in the uid are identied. In these cells, the velocity will be imposed using (2.56). Then, for each interface cell, the tangential velocity U 2 is evaluated in a point located on the same local normal, at a distance δ, as shown in gure 2.10. The position of point 2 is:

x 2 = x 1 + [δ -ϕ(x 1 )] ∇ϕ(x 1 ) ∇ϕ(x 1 ) . (2.57)
Figure 2.10: Implementation of the wall law. Found in [START_REF] De Tullio | Development of an immersed boundary method for the solution of the preconditioned NavierStokes equations[END_REF].

As a remind, ϕ represents the level-set function.Using U 2 , it is possible to obtain a rst approximation of the wall shear stress:

τ w µ U 2 δ . (2.58)
Therefore, we can evaluate the friction velocity:

u τ = τ w ρ νU 2 δ , (2.59)
by which the distance in wall units of the second point can be calculated and used to update u τ :

y + 2 = δu τ ν ⇒ u + 2 = f (y + 2 ) ⇒ u τ = U 2 u + 2 . (2.60)
This iterative process stops when convergence is reached. Once the nal value of the friction velocity has been calculated, the distance of the interface point in wall units is:

y + 1 = ϕ(y 1 ) ν .
(2.61)

The non-dimensional velocity is evaluated using the wall law:

u + 1 = f w (y + 1 ). (2.62)
Thus, the tangential velocity to impose at the interface cell is:

U 1 = u τ u + 1 . (2.63)
The choice of the distance δ is critical for the correct behaviour of the wall correction. Indeed, if the distance is too large, due to the parallelization (the choice has been done to have access to the rst layer of neighbours for each cell), the point where U 2 is evaluated may be outside the ghost layer. On the other hand, if δ is too small, U 2 would be calculated too close to the interface which is not accurate. So, we choose to impose δ equal to twice the biggest value of the level-set function at all the interface points. Finally, since the point x 2 is not necessarily at the cell-centre (where the variables are stored), velocity values should be interpolated in order to obtain U 2 . The interpolation is done by looking for neighbour cells through the closest node of the cell which owns point 2. Then Radial Basis Functions are used for the interpolation with this stencil. The wall function adjusts the velocity eld after the penalization step.

Chapter 3

Validation of the Numerical Tool

As described in the previous chapter, a numerical tool for solving incompressible Navier-Stokes equations on octree grids has been developed for two and three dimensional cases. Each main operator has been independently tested and compared with an analytical solution to conrm the space accuracy and the consistency. This chapter is now dedicated to the validation of the whole solver. First of all, the space consistency and accuracy of the 2D and 3D solvers have been evaluated. Then, the tool has been widely tested with benchmarks of two-dimensional cases found in the literature. The 3D code is then considered and a stability issue is reported. We found that the origin of this issue is the semi-Lagrangian scheme. Another scheme has thus been implemented and the tool was then tested with a benchmark of three-dimensional cases found in the literature. Finally, a comparison has been done between the octree code and a high-order Cartesian code.

1 Consistency and Accuracy of the Navier-Stokes Solver on Octree Grids

The rst validation performed for the whole Navier-Stokes solver is the analysis of the consistency and accuracy in space.

Taylor-Green vortex

The Taylor-Green vortex is an unsteady ow of a decaying vortex, which is considered in the present work since it has an exact solution in two spatial dimensions. A domain Ω = [-π, π] × [-π, π] is used and the solution of the Taylor-Green vortex is given by:

u = cos x sin y F (t), v = -sin x cos y F (t), (3.1)
where F (t) = e -2νt , ν being the kinematic viscosity of the uid and t the time. The pressure eld p can be obtained by substituting the velocity solution in the momentum equation and is given by

p = ρ 4 ( cos 2x + cos 2y)F 2 (t). (3.2)
Dirichlet boundary conditions have been imposed as ghost cells just outside of the borders of the domain thanks to analytical expressions. The computed solution as well as the mesh used for the analysis are shown in gure 3.1. Indeed, in order to get the consistency and the order of convergence in space of the whole solver, the computed velocity eld has been compared with the exact solution and the results can be seen on table 3.1. The time step used for this simulation Table 3.1: Norm of the spatial error and orders of the Taylor-Green vortex for 2D quadtree grid was ∆t = 10 -4 and the nal time was t max = 0.1.

A second-order accuracy is obtained for the L ∞ norm and a few more for the L 2 (and L 1 ) norm. These results allow a rst validation of our two-dimensional Navier-Stokes solver before going further with a comparison with data from literature in the next section.

Three-dimensional solver

Concerning the numerical tool in 3D, no exact solution exists for the whole computation of the Navier-Stokes equations to the best of our knowledge. Therefore, a partial analysis has been performed. An analytical function f exact (x, y, z) = sin(x 2 + y 2 + z 2 ) is considered and the Poisson equation is solved using the developed methods explained in the previous chapter:

∇ 2 f c = ∇.s, (3.3)
with s(x, y, z) being the analytical gradient of f exact . f c is so computed. The next step consists in computing with the developed methods the gradient of f c and comparing its values with the analytical expression s(x, y, z). The consistency and the order of convergence in space of these steps, corresponding to the projection and correction steps of the 3D Navier-Stokes solver, can be seen on Thus, we can validate the consistency of these parts of the solver. A rst-order accuracy is obtained for the L ∞ norm and 1.5 for the L 2 norm. A further analysis is indeed necessary to validate the whole solver and a benchmark of test cases using data from literature will be used.

2 Validation with Data from the Literature for the Two-Dimensional Solver

In this section, the 2D solver is widely tested thanks to benchmarks found in the literature. It was indeed necessary to completely validate the methods developed for the two-dimensional case before extending these methods to a 3D framework. The computation of the aerodynamic forces exerted on a body is necessary for this work to validate the solver. The aerodynamic forces have the following expression:

F = ∂Ω body -pI + νρ ∇ • u + ∇ • u T n dS, (3.4)
I being the identity matrix. The chosen method of computation uses the Lagrangian markers that dene the geometry of the body. A loop is carried out over the markers and for each marker, the cell that owns the marker is detected. Then, neighbours through the closest node of the cell to the marker are identied and a bilinear interpolation is performed with this stencil. Around the body, a Cartesian grid conguration is used in order to have the best accuracy possible in the computation of the aerodynamic forces. The pressure value is so obtained at the marker as well as the velocity derivatives, thanks to the bilinear expression whose gradient is straight forward to obtain. The simulations performed in this section used the PlaFRIM cluster, whose nodes are composed by a couple of 12-core Intel Haswell CPUs with a clock speed of 2.5 GHz and a shared memory of 128 GB (see https://www.plafrim.fr/).

Fixed Cylinder Case

The 2D cylinder is a test-case that is fully studied in the literature. It was chosen in the present work to focus on three dierent Reynolds numbers that correspond to three dierent ow regimes. Laminar ow with periodic vortex shedding is studied at Re = 100, transitionin-wake regime is observed at Re = 200 when instabilities lead to the beginning of formation of streamwise vortex structure and nally unsteady ow with low turbulence is seen at Re = 600.

Geometry and Mesh

All simulations performed in this section used a cylinder with diameter equal to 1. The size of the domain is 30 and the lateral boundaries are located at 15 diameters from the centre of the cylinder. The inlet is located at 7 diameters in front of the cylinder and the outlet is thus at 23 diameter behind. The mesh used can be seen in gure 3.2. Five dierent levels of renement are employed in this grid conguration. The aim is to simulate accurately the cylinder's wake. In order to choose the right compromise between accuracy and computing time, a grid convergence study has been performed.

Grid convergence study

In order to analyse the grid convergence, the case of a ow past an impulsively started cylinder has been chosen. Indeed, there exists a numerical solution of reference valid for short times as explained by Bar-Lev and Yang [5] and, especially, Ploumhans et al [START_REF] Ploumhans | Vortex methods for high-resolution simulations of viscous ow past blu bodies of general geometry[END_REF] studied the case of Re = 550. Thus, this ow has been studied in our work with three dierent meshes similar to gure 3.2 with dierent minimum level of renement. The characteristics of each grid and the computational time spent to simulate 6 s of ow are explained in table 3.3. The chosen time step for the simulations is adaptive and is computed at each time step. It depends on the size of each cell, on the velocity eld in that cell and on the CFL (Courant-Friedrichs-Lewy) parameter. This parameter has been chosen equal to 0.3.

The drag coecient history has then been compared with the analytical solution and can be seen in gure 3.3. This comparison allows to validate the developed solver. With a minimum tree level of 9 (which corresponds to around 270 cells on the diameter of the cylinder) the results obtained are very close to the analytical ones but this mesh implies a high computational time.

Therefore the grid with a minimum tree level of 8 has been used for the results described below.

The Strouhal number is a dimensionless variable describing oscillating ow mechanisms. It can be dened as:

S t = f D U ∞ (3.5)
with f the frequency of the vortex shedding, D the diameter of the cylinder and U ∞ the inow velocity at inlet. for the ow around a circular cylinder

Minimum

In order to characterize the transition-in-wake regime, the ow at Re = 200 has been simulated. Indeed, as explained by Williamson [92] this ow conguration corresponds to the beginning of formation of streamwise vortex structure. The ow can be seen in gure 3.6 with the isobars and the contour lines of z-component of vorticity and a von Karman vortex street can be observed.

The drag and lift coecient history of this simulation can be seen in gure 3.7. This gure was used to compute the averaged drag coecient and the Strouhal number of this ow. The results and the comparison with the other data from literature can be observed on table 3.5. The mean drag coecient shows a good agreement with the literature. As this two-dimensional simulation shows good results in simulation of a ow that could a priori be three-dimensional, that implies that these eects should be small enough, so that our two-dimensional solver is validated.

Flow at Re = 600

This ow regime at Re = 600 is faintly turbulent with appearance of three-dimensional effects (as explained by Williamson [91] [92]). The gure 3.8 presents the pressure eld and the z-component of the vorticity at t = 100. It can be observed that the vortices constituting the von Karman vortex street have a higher frequency than the previously studied ows. That

Airfoil Case

The solver developed has then been tested with the case of the ow around a 2D airfoil, that shows a more complex geometry than the cylinder. Data from literature have been used to compare with our solver as can be observed on table 3 Thus, the two-dimensional quadtree code has been validated for dierent ow congurations and for dierent geometries. The three-dimensional solver on octree grids is now studied and compared with data from literature for validation.

3 Stability Issue with the Semi-Lagrangian Scheme for Three-Dimensional Solver

While performing a low Reynolds simulation of the ow past a sphere on a hemispherically rened mesh, some numerical instabilities were detected in the proximity of the spherical leveljump, making the simulation diverge in a few seconds. An investigation has been performed to understand which part of the solver could imply such instabilities. Moreover, every operator was validated with an analytical expression. Every part of the Navier-Stokes solver has been tested and stabilized operators were implemented in order to stabilize the treatment of level jumps. But nothing changed and nally it turned out that only the advection term of the discretized Navier-Stokes equations was involved. The use of the semi-Lagrangian scheme as described in the previous chapter was unstable with some octree grids (not all kind of them). This instability is due to the fact that in the rst hand, we are using non-uniform grids and in the second hand, when we interpolate the velocity eld at the root of the characteristic curve the stencil used is not necessarily upwind. It can be noticed that, when stable, the semi-Lagrangian shows accurate results when compared to the data from literature.

An explicit Eulerian scheme has been implemented with a rst order upwind reconstruction. This was indeed the easiest way to proceed considering the remaining time of the thesis. The numerical tool developed constitutes actually a basic structure that will lead to further improvements. The prediction step (2.7) is so rewritten with the Eulerian scheme:

u * -u n ∆t + ∇ • (u n ⊗ u n ) = ν∆u * . (3.6)
A nite volume discretization has been used to compute the convective term on each cell volume

Ω c . ∇ • (u n ⊗ u n ) = 1 Ω c ∂Ωc (u n ⊗ u n ).ndS, (3.7)
where n represents the outward normal of the cell. Hence, we can write for each cell i of the domain:

∇ • (u n ⊗ u n ) i = 1 Ω c j (u n ⊗ u n ) j .n j ∆S j . (3.8)
The ux on the cell faces are then evaluated through a rst order upwind reconstruction:

(u n ⊗ u n ) j .n j = u i max(u j .n j , 0) + u e min(u j .n j , 0), (3.9)
where e is the cell pointed by the outward normal vector of the i-th cell, and the velocity at the face is basically:

u j = u i + u e 2 . ( 3.10)
Once this method has been implemented, the instability issue was totally solved as can be seen in gure 3.10. The three-dimensional solver has then been fully tested by comparing with benchmark from literature.

4 Validation with Data from Literature for Three-Dimensional Case

This section explains the work done in order to validate the developed 3D solver on octree grids. Two popular benchmarks have been used with the ow past a sphere at Re = 500 and the ow past a cylinder at Re = 3900. The ow past a sphere at low Reynolds number has been rstly studied to validate the numerical methods on a simple test case before dealing with higher Reynolds number for the case of the ow past a cylinder.

The computation of the aerodynamic forces that exerted on a 3D body has been implemented in the solver. The method is based on the Lagrangian markers (that dene the geometry of the body), but the method is slightly dierent than the two-dimensional case where bilinear interpolations were employed. A loop is carried out on the markers and for each marker, the cell that owns the marker is detected. Then, neighbours through the closest node of the cell to the marker are identied. A RBF interpolation is then performed with this stencil to directly get the pressure value at the Lagrangian marker position. The velocity gradients are also needed at these positions. The same stencil as previously explained is used. A Moving Least Square interpolation has been used as described in [START_REF] Gossier | Moving Least-Squares: A numerical dierentiation method for irregularly spaced calculation points[END_REF]. That method comes down to a localized least-squares procedure, in which the error E is expressed as:

E = N n=1 L n [u n -q(x, t)] 2 (3.11)
with N representing the number of neighbours identied, u n is the function to interpolate (the x-velocity for instance) where values are known at neighbour cell-centres. q(x, t) is the local rst-order polynomial interpolation whose components are the position dierences between each neighbour and the Lagrangian marker. The coecients of the local polynomial are thus directly the local gradients. L n is named the "localization parameter" and weights the contribution of the dierent points to the error. In the present work, we took L n = 1 Dist 3 equals to the inverse of the distance with a power of 3 after a preliminary study on this parameter. Minimization of E with respect to each of the coecient (local gradients) yields a 3x3 system of linear equations. The reader is refered to the report in [START_REF] Gossier | Moving Least-Squares: A numerical dierentiation method for irregularly spaced calculation points[END_REF] for the detailed linear system. Finally, solving these systems allow to obtain directly the gradients at each of the Lagrangian markers.

Sphere Case

This benchmark has been chosen to validate the operators developed in our solver with a low Reynolds ow conguration. All the simulations performed used a sphere with a diameter equal to 1. The size of the computational domain is 24 diameters and the lateral boundaries are located at 12 diameters from the centre of the sphere. The inlet is located at 8 diameters in front of the sphere and the outlet is thus at 16 diameters behind. An example of the mesh used can be seen in gure 3.11. Four dierent levels of renements are employed and a grid convergence study has been performed with dierent minimum level of octree. The characteristics of each grid can be seen in table 3.8 as well as the conditions of the simulations. The PlaFRIM cluster has been employed for these simulations. An adaptive time step has been used.

Minimum tree level

Minimum size The choice of the CFL number is very important and the results show really dierent behaviours depending on its value. The gure 3.12 describes the time history of drag coecient with the 3 dierent octree grids and with a CFL condition of 0.5 compared to a CFL of 0.3. The higher CFL implies drag coecients too high compared to existing value in the literature. Therefore for the future simulations a CFL condition of 0.3 has been imposed.

Number of cells on the diameter

Global number of cells

Number of processors

The averaged drag coecient obtained is in line with results found in literature as remarked in table 3.9. With the minimum level of 7, it seems like the grid convergence is reached since the

Re

Authors

Cylinder Case

In the following sections, the simulations have been performed with the supercomputer Occigen 1 . The ow past a cylinder benchmark at Re = 3900 has been used to test the high Reynolds behaviour of the octree Navier-Stokes solver, coupled with the Wall Modelled Large Eddy Simulation approximation (WMLES). Since the boundary layer is laminar at the considered Reynolds number, the Reichardt wall function is not suited. Therefore, a linear correction is implemented: u + = y + , which can be considered as a sort of laminar wall function. As a reminder, u + represents the non-dimensional tangential velocity and y + is the dimensionless wall distance. This case constitutes a rst test and validation of all the models developed for the 3D code before studying a higher Reynolds ow in the next section. The ow regime at Re = 3900 is still challenging since transition waves appear in the free shear layer in the wake. One diculty with our solver to simulate properly this test case is the fact that we are constrained to use cubic domains. Thus, if we want to have a large wake visualization, that implies also large lateral borders. As we are studying an innite cylinder, no coarsening can be done in the direction of the cylinder length and a high number of cells arrives quickly. A compromise should be accepted between the accuracy in the wake of the cylinder and the number of cells. The chosen grid has the following characteristics: Several experiments and computations have been carried out for the ow past a circular cylinder at Re = 3900 which will allow to compare the dierent data from literature with the developed solver. But we should be aware that this numerical tool, allows to have an interesting compromise between accuracy and computational time but we do not expect to be as accurate as high-order numerical methods. A simulation of 90 s of ow has been performed with 2880 processors. The associate computational time was around 72 hours. An adaptive time step has been used with a CFL condition of 0.3. For the rst comparison, we used the computation of the averaged drag coecient as can be seen in table 3.10. Our solver shows to be in line with data from literature for the mean drag coecient. The pressure coecient on the cylinder has then been studied and compared with experimental values of Norberg [START_REF] Norberg | Eects of Reynolds number and a low-intensity free-stream turbulence on the ow around circular cylinder[END_REF]. This computation has been done after a preliminary simulation for removing the initial transient. Pressure values have been computed at the middle of the cylinder and are averaged over 9 vortex sheddings (corresponding to around 40 seconds of simulations). We were indeed limited by the computational time of the simulation since we had a limited number of hours of simulation on the supercomputer Occigen. The results obtained can be seen In the same way as for the pressure coecient, the numerical results are obtained after a preliminary simulation and by performing an average over a period of 9 vortex sheddings. Figure 3.15 shows the wake prole at three dierent locations (x = 1.06, x = 1.54 and x = 2.02). The results from the octree solver are in-line with the experimental data of Parnaudeau et al. [START_REF] Parnaudeau | Experimental and numerical studies of the ow over a circular cylinder at Reynolds number 3900[END_REF] for the position close to the cylinder (at x = 1.06). The results obtained for the two other positions, farer from the cylinder, are then closer to the experimental data from Lourenco and Shih [START_REF] Lourenco | Characteristics of the plane turbulent near wake of a circular cylinder, a particle image velocimetry study[END_REF]. Indeed, a level jump is present in the wake at x = 2 which should imply a loss of accuracy in this area. Also, the period used for averaging the wake prole data could not have been very wide due to computational limits.

Re

Authors

Finally, the results obtained with the developed solver are in accordance with the data from literature. They are not as accurate when looking at the velocity prole in the wake but this is probably due to the numerical dissipation implied by the rst-order accuracy of the numerical tool. The next section will focus on a higher Reynolds number ow and take particular attention to compare the simulations performed and the results obtained with a high-order solver on Cartesian grids.

Comparison between Cartesian and Octree Solvers

The ow past a cylinder benchmark at Re = 140000 has been used to compare simulations performed with the developed octree solver and a Cartesian code that was developed before the present work. This solver has been developed and validated in-house by the MEMPHIS Team of INRIA and is called NaSCAR. High-order nite dierence model was implemented with an immersed boundary method. Considering the simulation with the octree code, the same domain size and mesh conguration as previously described for the case at Re = 3900 has been used. The Cartesian code is not submitted to the constraint of having a cubic computational domain. Therefore the size of the domain chosen for the simulation was the following: The grid used has 2400 × 1600 × 200 cells respectively in the x, y and z-directions which implies a total of 768 millions of cells. Around 100 cells are located on the cylinder diameter which is a few more than with the octree grid conguration. If we resume, the octree conguration has a domain around 12 times bigger than the Cartesian domain, the number of cells is around 5 times lower than the Cartesian grid and the number of cells on the cylinder diameter is 1.18 lower than for the Cartesian conguration. The number of unknowns is not the only factor of computational time since the performance of the solver has also a signicant impact. The Z-ordering used to deal with the discretization on the octree grid implies indeed matrices less well conditioned than the discretization with the Cartesian solver NaSCAR. No test has been performed with exactly the same domain, same grid, same ow and same computational machine but we assume that the Cartesian solver is around 3 times faster to solve the Navier-Stokes equations than the octree solver. For this case, for example, the machines used for the simulations were not the same and especially the CPUs had dierent clock frequencies. The simulations on the Cartesian mesh with NaSCAR have been performed with the supercomputer Turing2 .

The rst comparison performed is visual and can be seen in gure 3.16. It can be seen that small structures of turbulence are captured with the high-order Cartesian code, whereas the octree code can't resolve that. But the ow patterns seem however to be correct with the octree code. Further comparisons have been performed by comparing both codes with data from literature. The averaged drag coecient has been computed for both simulations and compared with data from literature. This can be observed on table 3.11. If looking at the experimental value and at the most recent simulation performed with LES done by [START_REF] Kim | Large eddy simulation of turbulent ow past a circular cylinder in subcritical regime[END_REF] [START_REF] Kim | Large eddy simulation of turbulent ow past a circular cylinder in subcritical regime[END_REF], the result obtained with the octree code is slightly more accurate than with the Cartesian code. Both results are in line with the data from literature.

Other comparisons have been performed to evaluate the accuracy of the solvers in the wake of the cylinder. For this purpose, experimental results from Cantwell (1983) [START_REF] Cantwell | An experimental study on entrainment and transport in the turbulent near wake of a circular cylinder[END_REF] have been used. This has been computed after preliminary simulations in both the Cartesian and the octree codes to remove the initial transient. The following results have been averaged over 5 vortex shedding (corresponding to around 20 seconds of simulations) for both the Cartesian and the octree code. Indeed the computational time needed was very high for the Cartesian case, therefore it has been chosen to perform the comparison on the same number of vortex shedding. The results can be observed in gure 3.17. At the location x = 1 behind the cylinder, the results from Cartesian code are closer to the experimental values. However, farer from the cylinder, the results from the octree code seems to better represent the experimental work. It was not possible to nd the number of vortex shedding used to perform the mean experimental wake proles but we expect that a very large time has been used. We know, that averaging on 5 vortex shedding is very little. Therefore, for the octree solver, other simulations have been performed to see if the wake prole remains stable. The results obtained after 5 other vortex shedding can be seen in gure 3.18. It can be observed that the streamwise velocity at x = 1 does not change much. However, at x = 3, the simulated results are farer from the experimental one, while staying in-line, especially compared to the results seen previously from the Cartesian solver. The normal velocity shows a signicant dierence between the two dierent averages. We don't know yet how to explain this behaviour, thus, an other simulation has been launched with the Cartesian solver with another 5 vortex shedding to see if the same behaviour occurs.

In the mean time, a very long simulation is performed with the octree code in order to average with a high number of vortex shedding. Perhaps this eect is due to a long period oscillations and necessitates a long time averaging to be compensated.

Finally, the code developed allows to get an interesting compromise between computational time and accuracy. The benchmark of test cases well described in literature allows to show that with the octree solver accurate results can be obtained when dealing with a reasonable number of degrees of freedom. The cylinder case has involved all the numerical models developed. The constraint of having a cubic domain prevents us to have a grid ne enough in the wake. Thus, the results obtained are not as accurate as a high-order solver, especially the Cartesian code developed in-house that is able to catch the small structures of turbulence. But this chapter allowed to validate the developed numerical tool with non-trivial test cases at high Reynolds number ows. The implementation of the application of our solver to the case of a rotating wind turbine will be described in the nal section. measured on a wind blade as well as the measurements of upstream wind is reported. A description of the technical solutions allowing to achieve this, is done. Then, the on-site sensor's installation process is explained and nally, the analysis of the data is described.

1 Specications

Global specications for the experimental work are rst to use a Jeumont J48 wind turbine because this is the only model on which we have access to geometrical and structural data. Some characteristics of this wind turbine are a hub height of 46 m; a rotor diameter of 48 m; and blades measuring 23.5 m long. Then, the duration of these measures should allow obtaining enough wind data, which should include several strong wind conditions. So, the planned period of measurements will be between 4 and 7 months. Indeed, when the specications were written, we had no idea of the cost and of the resistance along time of future sensors on the wind blade.

In order to reach the previously explained objectives of this experimental work, discussions were held with INRIA to enable identifying the data of interest to measure.

Wind Measurements

The measurement of the wind was considered as the most important data to be obtained from the very start of the work. Indeed, this data will be used as a key input for the numerical simulations. This is needed to represent the wind as it arrives at the wind turbine. Figure 4.1 shows an example of computational domain for numerical simulations. Thus, the measurement should be done in a wind sector with a high occurrence in the wind rose of the site. In order to have representative data of the wind arriving at the turbine, no perturbations from the terrain should occur between the measuring tool and the wind turbine. So, the terrain should be quite at and without many trees and obstacles. Installing the measuring tool in the wake of another wind turbine should also be avoided.

Another requirement is to enable estimation of the vertical wind prole. Indeed, the wind velocity prole in the vicinity of the ground generally shows a strong gradient (in the direction normal to the ground). This phenomenon is known as the atmospheric boundary layer and depends on many parameters, such as the roughness of the surface, the weather conditions (pressure and temperature) and also thermal stability changes between day and night. So, several measurement points are needed at dierent heights. In order for the future numerical simulations to be accurate and similar to the real experiment, the turbulent ow should be measured.

Specications

Moreover, the wind during a gust phenomenon is characterized by high variations in wind velocities in a few seconds as explained in the IEC standard [1]. Thus, the acquisition frequency of wind data should be small enough to accurately characterise these changes.

Blade Data

The main goal of the wind turbine data collection is to observe the aero-elastic behaviour of wind blades. The decision was made to collect data from the maximum number of measurement points, given the budget in the AEROGUST project. Thus, it was decided to fully instrument one wind turbine blade (instead of all three blades).

A global requirement on instrumentation of the wind blade is, rst, to minimize the impact of measurement tools on the wind blade. Indeed, the wind turbine is supposed to keep working well after the experimental phase. So, for example, no hole is allowed to be made. Then, the Valorem company does not have the resources to develop an inquiry system to get data back from sensors to the oce. Thus, we looked for a complete technological solution with regards sensors supply, installation, maintenance and communication system via a GSM modem to collect the data. In order to have experimental data showing the aero-elastic behaviour of a wind blade, the measurement of the blade deformations will allow to check the structural beam model of the blade and to quantify its structural behaviour. To observe the aerodynamic load on the wind blade, the measurement of pressure of air on the blade will be of signicant interest.

The aim of the blade deformation data will be, rst, to check that the beam model used in the numerical simulations is in accordance with the real deformations of the blade and second, to observe the eects of the wind on the blade deformations. Strain gauges are considered for this measurement of blade deformations. For the purpose of validating the beam model of the blade, a requirement is to put several strain gauges on the blade root, since the strain at the root is the sum of all strains along the blade. The measurement of experimental pressure data will be used for comparison between experimental data and numerical simulations.

Moreover, we realized that the correlation of the wind on the blade and the analysis of the blade deformation will need to be spatially synchronised. Indeed, as this will be the rst access of such structural data, the dierent kinds of deformation that the blade undergoes in reality is to date relatively unknown, while being related to its position on the rotation circle of the wind turbine (for example, when the blade passes along the tower, the behaviour of the deformation is probably not the same as when the blade is clear of the tower). So, another requirement for the experimental work is to sense the position of the blade relative to the tower (rotational position).

The operating range encountered by the sensors necessitated writing specications to our suppliers. For the blade deformation, data from the blade manufacturer was used to determine the expected range of deformation. Thus, the local deformations of the blade are of the order of -1.5 % to 1.5 % and the tip deection of the blade can reach 2 m. For the pressure, the Bernoulli equation was used and applied to the approximation of the mean air density and expected normal wind velocity range. So, the operating range of the pressure sensors is of the order of 1/10 th of the atmospheric pressure (on average 104 P a). The positioning of the sensors was not decided during the specication step, since it was unclear how many measurement points the budget would allow.

A requirement was written to specify the targeted accuracy of experimental data. Indeed, the numerical results from our computational code are expected to have an error of the order of 5 %, so a similar accuracy is required from the experimental results. The nal requirement for all the blade data is to have the same acquisition frequency as the wind data in order to ease the synchronisation and the correlation between both kinds of experimental data. Now, the manner in which the technical solutions could meet the specications will be explained in the next section.

2 Technical Solutions for Measurements

Meteorological Mast

The aim is to have the most accurate measures and to be in accordance with the budget.

Valorem has experimented 3 dierent ways to measure wind velocity and turbulence at dierent heights. The measurement with a meteorological mast, with a SODAR (SOnic Detection And Ranging) which measures the scattering of sound waves or with a LIDAR (LIght Detection And Ranging) which measures the reected pulses of a laser light with a sensor. Valorem owns several meteorological masts and Sodars and the cost of a Lidar is very high. The IEC 61400-12-1 standard [1] is the reference for wind assessment studies. The standard prescribes the use of a meteorological mast in such studies. Moreover, the Sodar has a maximum acquisition frequency of several seconds which is not adequate for our specications. That's why the nal choice was the use of a meteorological mast (met mast).

Valorem had a 50m height meteorological mast available. As the height of the hub of the wind turbine is 46 m, an extrapolation of the wind to the top of the blades will be done thanks to the measured wind prole. All the sensors will be oriented to minimize the wake eect from the mast according to the IEC 61400-12-1 standard. According to this same standard the met mast should be located between 2D and 4D (D being the turbine rotor diameter) upstream from the wind turbine, ideally at 2.5D, which corresponds to approximately 120m.

Then, work has been done to dene the setting-up of the met mast. Indeed, requirements are to (a) nd a at terrain without perturbations (trees, wakes of other wind turbines. . . ) between the wind turbine and the mast in the main direction of wind; which (b) the owner of that terrain accepts grant for this purpose. A wind rose denes the main wind directions. It has been obtained thanks to a wind measurement campaign by Valorem in another terrain located a few kilometres away of our site. As can be seen in gure 4.2, the wind plant is located on a hill so the met mast should also be on the same height. Indeed, putting the met mast upstream from wind turbines n°4 and 8 should be avoided because of the wake of other turbines. Turbines n°1, 2 and 7 have strong terrain variations so are not good choices either. This leaves wind turbines 3, 5 and 6 as possible options.

Going forward, a site-inspection was conducted to conrm the possibility for the installation of the met mast and to investigate whether the owners of elds would agree. Compensation was planned for owners and farmers of the terrain where the met mast would be placed. But as wind turbines are not always well accepted in France, it was not easy to obtain an agreement from the owners to install the met mast. Finally, a terrain was found upstream from wind turbine n°5. The next step was then to obtain an administrative permit for the setting up of the met mast as will be explained in the next section

Blade Measurements

The specications for the blade measures have been sent to several companies. The only viable solution was suggested by a company named Febus Optics, based in Pau (France). They suggest the use of optical bre with Fiber Bragg Gratings (FBG) sensors. Achieved by creating a periodic variation in the refractive index of the bre core, FBG are short segments of an optical bre that reect particular wavelengths of light and transmits all others as can be seen in gure 4.3. This wavelength can be measured by injecting an incident light having a broad optical spectrum and analysing the wavelength of the reected light with an OSA (Optical Spectrum Analyser). It can be used as a sensor because any change in the period of the refractive index will induce a change in the specic wavelength (so-called Bragg wavelength) reected by the FBG. Thus, strain or temperature will induce a measurable wavelength change. Many types of sensors can be addressed by this technology as soon as the body of the sensor converts the measure in wavelength change in the FBG.

Moreover, using dierent Bragg gratings with dierent wavelengths, it is possible to connect several sensors inline. That means that every sensor will be connected to the others in the bre and the dierent bres will then be connected to an inquiry system located inside the rotor. More precisely, on the rotating part of the rotor, so that no problem occurs with the bre arriving in the hub (it would be very complicated to pass from the rotating blade to the xed part of the rotor). The company also suggested the installation of an inclinometer inside the inquiry system to measure the angular position of the instrumented blade. This technological solution is very non-intrusive since all the sensors and the bres would be stuck to the blade surface. The bre is also exible which allows a wide range of sensor locations. The inquiry system will collect the data and send them via 3G to the oce. Finally, this technological solution ts all the specications and the cost is in accordance with the budget. This was therefore the technological solution of choice. Febus Optics works with another French company, based in Brittany (France) close to the wind plant, which manufactures the optical bres, pressure sensors and strain gauges.

Position of Sensors on the Wind Blade

Concerning the number of sensors (for pressure and strain) to put on the blade, the main constraint was the budget. Indeed, the more sensors installed, the better our experimental work. Strain gauges were less expensive than pressure sensors, so more of these were chosen. This allowed the purchase of 10 strain gauges and 6 pressure sensors.

For strain gauges, as explained previously, one of the goals is to check the beam model. So, strain gauges should be put on the blade root, where the strain is the sum of all the strains on the blade. It was decided to place 4 strain gauges every 90°at the cylindrical blade root. Then, Figure 4.5: Position of the strain gauges on the blade to optimize the use of the remaining strain gauges, their positions were chosen to be where only longitudinal strains occurs. Indeed, if we would like to measure the strains in several directions, several strain gauges would be necessary at each location. The other locations have been chosen to have a kind of equipartition along the length of the blade. The locations of the strain gauges can be seen on gure 4.5. At section n°2 we put the gauges at the same locations as the blade root i.e. at the leading edge, trailing edge, and on the sparcaps at the pressure and suction sides. At section n°3, two strain gauges were located on the sparcaps at the pressure and suction sides.

As we only had 6 available pressure sensors, their locations were chosen in order to obtain a pressure prole around a section, and in a manner similar to experiments conducted by the NREL (National Renewable Energy Laboratory) on a small wind turbine [START_REF] Simms | -IV: Test congurations and available data campaigns[END_REF]. The sections where the pressure sensors would be placed was chosen as a compromise between two requirements. First, the sensors should be located far from the nacelle since the nacelle geometry is not well known and can't be reproduced accurately in our simulation. Then, we cannot put the sensors where the chord is too small, otherwise the geometric perturbations induced by the sensor on the airfoil would become signicant (the sensor has a thickness of 1 cm, a length of 5 cm and a width of 2 cm). The chosen location for the pressure sensors can be seen on gure 4.6.

In the next section, work done before and during the installation on site of the met mast and of the sensors on the blade will be reported.

3 On-site Installation

Meteorological Mast Installation

Prior to the met mast installation, some administrative tasks should be performed. First, the civil aviation services must be consulted and should agree. Then, an administrative permit (named "preliminary declaration of work") should also be delivered by the city council. This le should include geographic location of the future met mast, photographs of the mast in its environment, drawings of the mast with its dimensions (2D and 3D) and how long the mast The nal step is to submit a DICT (Declaration of the intention to begin work). It consists of checking than the future mast will not impact on power lines, gas pipes, water network, etc. All this was subsequently done by Valorem for the project.

The nal position of the met mast can be observed on gure 4.7 and it can be observed thanks to the wind rose that the wind turbine and the met mast are aligned with the main wind direction occurrence on the site.

A brief description of the installation process for the met mast follows:

The mast is composed of several tubes that t together on the ground. 8 anchoring points are made in 4 orthogonal directions (in order to have a high stability of the mast).

Installation of wind vanes, cup anemometers, logger, photovoltaic panels, nocturnal runway lights is done on the ground.

The mast is raised thanks to a motorised winch.

Adjustments are made on the guide wires so that the mast stands vertical and straight.

The certication of the cup anemometers is done thanks to wind tunnel measurements. The calibration laboratory certied that the value of the measured wind lies within the assigned range of values with a probability of 95%.

Sensors on Blade

For on-site installation of the blade sensors two choices exist: using steeplejacks or an aerial work platform. Although an aerial work platform has a signicantly higher cost, the time needed with steeplejacks for installation would be too long and so too complicated to implement the sensors on the blade.

The supplier of the sensors had no previous experience in wind blade instrumentation. Thus, preliminary checks were done. Indeed, Valorem has available nearby its oces, a 48m long wind blade of the same model than the blade that will be instrumented. The blade is located on the ground. A test of implementation of sensors on the blade was therefore done at Valorem before the real aerial installation. This test revealed that the use of the high adhesive tapes as rstly planed wasn't a good option. Indeed, thanks to the previous experience from Valorem, the fact that a step exists between the adhesive tape and the blade, renders the surface not smooth enough, and is likely to cause a problem. The blade is submitted to signicant load and apparent wind, so, these kinds of steps with an obstacle on the blade could lead to the loosening of the tape

The methodology suggested by the supplier was then the following:

Sanding of the painting of the blade where the sensors and the optical bres will be located from the top (the nacelle) to the bottom along the blade.

Cleaning of this surface with acetone.

Passing the bre through the hub of the wind turbine for its connection to the inquiry system.

Positioning of the bre and/or sensor and gluing with epoxy glue.

Protecting with low adhesion wove tape.

Painting with protecting polyurethane paint designed for wind blades. Once the met mast has been installed and it has been checked that the wind data was correctly collected, the installation of the sensors on the blade began. Necessary conditions for installation are complete absence of both rain and strong wind. Indeed, rain is a problem for gluing with epoxy and a strong wind prevents the use of the aerial work platform. The prevailing weather was therefore monitored prior beginning the installation. We had to be sure that as soon as we began the implementation, we would be able to nish it. Finally, a half day delay occurred because of strong wind and the whole installation lasted almost 3 days. A view of the work on the aerial platform can be seen on gure 4.10. In gure 4.11, the 2 lines of sensors going along the pressure side and the leading edge of the wind blade are shown. The setup of the 4 cables is:

Pressure side: 3 strain gauges, 1 pressure sensor.

Leading edge: 2 strain gauges, 2 pressure sensors.

Suction side: 3 strain gauges, 1 pressure sensor.

Trailing edge: 2 strain gauges, 2 pressure sensors.

As a reminder, the box containing the inquiry system is located on the rotating part of the rotor and contains an inclinometer to measure the angle of the blade. Finally, organization of the on-site installations was successfully conducted for both the met mast and the sensors on the blade. Work has then be done with the experimental data in order to identify dierent gust conditions in the eld and to analyse the eects on the blade deformations. The measurement campaign lasted 9 months from April 2017 to January 2018. At the end, a huge number of experimental data have been collected. Considering the wind speed measurements, 26 939 646 data points exist for each 5 anemometers. For the blade data, 11 806 822 measurements are available for each of the 16 sensors (10 deformation data and 6 pressure data). This corresponds to an availability for the meteorological mast of 88 % and 43 % for the blade's sensors. Indeed, several issues arose along the campaign that implied some data loss. The methodology used will so be described and then the most important outcomes of the study will be explained.

Methodology

Considering the amount of data, the rst step in the analysis of this data was ltering to ease the handling of data. Indeed, the wind speed measurement should be representative of the wind arriving on the turbine so the wind direction should be lined up with the met mast and the wind turbine. Direction of the wind has so been ltered at ± 10°around the optimal direction. We also wanted to have a rotational speed of the rotor representing a normal behaviour of the wind turbine. That means for example to avoid the start of the turbine. Therefore, the blade rotational speed has been ltered to keep only the value between 11 and 25 rounds per minute (25 rpm corresponding to the maximum rotational speed of the turbine).

Once, this work of ltering was performed, the rst idea was to gather all the remaining data along time and to see if a correlation between wind speed and sensor's measurement could be Figure 4.12: Longitudinal deformations from sensors at dierent locations as a function of wind speed observed. Figure 4.12 shows the deformation data from the strain gauges located at pressure and suction side of the blade as a function of wind speed. The inuence of the wind speed intensity in the values of deformations can be seen but no signicant correlation has been highlighted with this methodology. Therefore, we change our strategy for the analysis of experimental results.

Another idea was to nd a methodology to identify the highest deformations occurring on the blade and look then how it could be related to the wind speed. Indeed, we are nally interested in knowing what produces the highest loads that exerted on a wind blade. However we didn't nd an adequate strategy to do that, since the strain gauges data show high variations along time whatever the wind is. Also, at Valorem, we are not specialists enough to characterize the deformations that would imply high loads on the blade. The choice was thus to focus on wind data rst and to identify gust events.

Identication of "gust events"

The IEC standard [1] denes a theoretical gust as a high variation with a Gaussian shape of wind speed on a period of several seconds (around 10 s). Therefore, the basic idea tested is to compute a moving average of wind speed on 10 s period. Then, on each period, the computation of the variation of wind speed is done (V max -V min ). A second method has also been tested, that consisted in using denitions of the "Extreme Operating Gust" (EOG) in the IEC standard. This concept denes an extreme wind condition due to a gust and an empirical equation exist to compute a value of the EOG from 10 minutes averaged wind speed data (and other parameter as the turbulence, the diameter of the rotor, ...). Using this equation, Extreme Operating Gusts have been computed from 10 minutes averaged wind speed data, which gave us a curve with a length of 10 seconds describing the EOG during each 10 minutes. A comparison of this theo-Category 1 Category 2

Number of gust events 28 10

Example Table 4.1: Characterization of gust events retical gust to the 1 second experimental wind speed data is then done to nd the correlation between the two. The two dierent methods have been tested on 1 month of experimental data and the method based on the IEC standard (i.e. based on the EOG denition) showed poor correlation coecients. Indeed, we analysed after, that this equation denes an extreme event but we had a few chance that such an extreme event occurred during our experimental work. We focused thus on the rst method (based on the variation of wind speed). It was noticed that, only doing V max -V min only nds gusts with high wind speeds, whereas normalizing with the minimum velocities V max -V min V min keeps high variations whatever the velocity is. These two variants of the method have been used, since we are interested in every kind of gust event.

This methodology has been applied to the whole period of available experimental data. The threshold used for our parameter to dene whether the variation corresponds to a gust or not, has been chosen so that, at the end, only a few events are left. Indeed we had no expertise to dene such a value otherwise. Finally we identied 38 gust events on 1 minute periods. It was noticed that the gusts can be characterized and classied into two dierent categories. We observed that two dierent behaviours of the gust event could occur, as summarized in table 4.1. On the one hand, we can see gusts that consist in a high increase of wind speed followed by a decrease to come back to the initial value of wind (category 1). This category shows gust events quite similar from what is described in the IEC standard. On the other hand, we can also observe a high increase of wind speed but then, no decrease occurs right after (category 2).

Comparison of deformation data between a gust event and with weak wind

A test has been done to compare the deformation data during a gust and with weak wind in order to validate that dierent behaviours occurred. First-of-all it was needed to reconstruct the time series of wind speed data measured at the meteorological mast so that we can evaluate when the gusts arrive on the wind turbine. This was done basically by transporting the wind speed along time. It implies of course uncertainties on the wind data but it gives an idea when the gust arrives at the turbine. The dierent behaviours of the deformation data during a gust or a period with weak wind can be seen in gure 4.13 with an example for one strain gauge located at the leading edge at 7 m from the blade root. When the wind is weak, it can be observed that the oscillations dening the deformations of the blade have a regular amplitude and frequency. The frequency has been computed with the time between two successive oscillations. However, during a gust, variation of the frequency of deformation can be seen around the gust for the sensors located at the leading edge of the blade. The sensors located at the leading and trailing edges show similar results, and, in a same way for the sensors located at the pressure side and suction side, which makes sense. The sensors located in the pressure and suction sides are also subjected to variation of the amplitude of data around the gust and less eected by variation of frequency. These are so rst remarks made from gures plotted for several events one by one.

Synthesis of the Experimental Work

In order to synthesize the results obtained from the experimental work and to better understand the eects of gusts on the wind blade, we came back to our very rst test that was to gather the data along time as a function of several parameters and see if a correlation of gusts and deformations could be seen. But, we only used data of the identied gust events and instead of the wind speed measured at the meteorological mast, we used the time series reconstruction. Therefore, as we noticed above that the frequency of the deformation seemed to increase when dealing with gusts, the mean frequency during periods of weak wind and during gusts have been compared. Indeed, we identied periods of weak wind (same number of events as gust events) to allow a serious comparison. The results can be seen in table 4.2. We can thus conrm that the frequency of the deformation is higher for every sensor in the case of gust events. Concerning the amplitude, the analysis is more dicult to be conducted. Indeed, the mean amplitudes are higher when submitted to weak wind compared to gust. An explanation could be that, as the frequency is higher during gust events, and the acquisition frequency of the data being 1 s, it is likely that several data points at high amplitude are missing in the data collection during gusts.

f
An interesting gure is the plot of the deformation as a function of the position of the blade. This can be seen in the gure 4.14 for the cases of gusts and during weak wind. The angle of 0 corresponds to the position where the blade is along the tower, the tip blade being at its closest position to the ground. The blade rotates clockwise. It can be observed that during regular and weak wind, the angle of the blade has a signicant inuence on the deformations located on the trailing and leading edge (same behaviour for both). This is also true at the pressure or suction side as we can see that points are less diused than during gusts but slightly less signicant. It can be remarked, that during gusts, the deformations of all the sensors seem to increase when the blade is going down. This phenomenon does not occur when the wind is weak. We can also notice, as expected, that for the sensors located at the blade root, the deformations at weak wind are very small.

Another remark is that, when plotting the deformations as a function of wind speed at the wind turbine (reconstructed from the measurements at the meteorological mast), no signicant inuence can be seen. The same kind of gure as the very rst test (see above gure 4.12) is obtained, whereas we are only dealing with gust events. This is an unexpected result since we imagined that the wind speed intensity would have an impact on the amplitude of deformations. Another hypothesis is that, since the position of the blade has a signicant impact on the deformation of the blade, we should subtract in a way this eect of the position (of the rotation of the blade, nally). Perhaps, this could lead to highlight the inuence of wind speed on deformation data. The results presented in this section constitute preliminary outcomes. Several gust events have been identied, that can be used as test cases for our numerical tool. Gusts that can occur on an operating wind farm have been characterized. This experimental work allowed a gain of experience for Valorem in terms of data treatment and knowledge on wind gusts. This analysis can be pursued in order to show the inuence of wind intensity on the deformations. Indeed, the inertial eects can be subtracted from the deformation data in order to highlight the eects due to wind.

Chapter 5

Towards Application to the Study of Wind Gust Eects on a Blade

The aim of the thesis is to simulate the ow around wind turbine blades, especially when submitted to gusts. As we have available several experimental data with deformations at dierent locations and pressure, it is relevant to implement a structural model in our numerical tool. This can so be used to calibrate our numerical tool and to compare the data with the simulations. Using a computational domain with the whole wind turbine and simulating the three blades, even with an octree grid, would require a huge amount a computational time. This would indeed imply a wide computational domain and ne cells around the three blades. Therefore, the approach chosen in the present work, is to simulate only one wind blade in a rotational framework. The rotational framework allows to use a signicantly smaller computational domain compared to the domain that would require the simulation of the rotation of the blade. This simple approach is a base but, in the future, the whole wind turbine could be simulated with the tower and the nacelle. This chapter describes the structural model used for the wind blade and the way it has been coupled with the Navier-Stokes solver. Then, the rotating framework and the required modications are explained. An adaptive mesh renement process has then been implemented and the results obtained are presented. Finally the simulations including all the numerical methods implemented are described with the results obtained.

1 Coupling with the Solid Model of the Wind Blade Apx48

Structural Model

Wind turbine blades are very slender structures with an internal conguration very similar to that of airplane wings, their structural response can therefore be described with a reasonable level of accuracy using beam models. In this work the model reported in [START_REF] De Frias Lopez | A 3D nite beam element for the modelling of composite wind turbine wings[END_REF] is used. The spatial discretization of the elasticity problem relies on the nite element method. The considered element possesses two nodes, each with six degrees of freedom: three displacements and three rotations.

Let X-Y-Z be the global reference system as can be seen in 5.1:

X axis corresponds to the beam axis, Y axis is parallel to the lagwise movement direction (apwise bending axis), Z axis is parallel to the apwise movement directions (lagwise bending axis).

Figure 5.1: Blade cross section with suggested reference system including denition of apwise and lagwise bending axes. Source [START_REF] De Frias Lopez | A 3D nite beam element for the modelling of composite wind turbine wings[END_REF].

The kinematics of the beam are described using the Euler-Bernoulli theory for the bending combined with the Saint-Venant hypothesis for the torsional motion. In the global system, the displacement eld is:

     u(x, y, z) = u t (x) + (z -z t )θ y (x) -(y -y t )θ z (x), v(x, y, z) = v c (x) -(z -z c )θ x (x), w(x, y, z)
= w c (x) + (y -y c )θ x (x), (5.1)
where y t and z t are the coordinates of the Tension center (T), y c and z c the coordinates of the Shear center (C), u t is the axial displacement of T, v c and w c are the bending displacement of C. Let X'-Y'-Z' be the reference system centred in T and parallel to the global one, the deformation eld is:

     xx = u t,x + (z -z t )θ y,x -(y -y t )θ z,x , xy = -(z -z c )θ x,x , xz = (y -y c )θ x,x . (5.2)
Since modern wind turbine blades are made of glass-bre and reinforced with lightweight carbon bres, the stress-strain relation is that of an anisotropic materials:

    σ xx σ xy σ xz     =     E 11 E 12 E 13 E 21 E 22 E 23 E 31 E 32 E 33         xx xy xz     .
(

The internal reactions are dened as:

                         N = A σ xx dA, M x = A [(y -y c )σ xz -(z -z c )σ xy ]dA, M y = A z σ xx dA, M z = A y σ xx dA. (5.4)
Substituting eqs. (5.2), (5.3) in (5.4) and integrating on the cross sectional area, the following system of equations is obtained:

       N M x M y M z        =        EA AT 0 0 AT GJ F T -LT 0 F T EI y y -EI y z 0 -LT -EI y z EI z z               u t,x θ x,x θ y,x θ z,x        . (5.5)
where:

EA is the axial stiness;

GJ is the torsional stiness;

EI y y , EI z z , EI y z are the bending stinesses, respectively around the Y' axis, Z' axis and the ap-lag coupling term;

AT , F T and LT are respectively the axial-torsion coupling, the ap-torsion coupling and the lag-torsion coupling, they are strongly dependent on the orientation of carbon / glassbres.

The balance equation of elasticity problem can be written using the virtual work principle:

W int = W ext , (5.6)
where W ext is the work of the external forces and W int that of the internal forces:

W int = L 0 (u t,x N + θ x,x M x + θ y,x M y + θ z,x M z )dx. (5.7)
The spatial variation of the displacements is described using the so-called shape functions:

u = N s , (5.8)
where {u} is the displacement vector:

u = u t v c w c θ x θ y θ z T , (5.9)
[N ] is the shape function matrix:

N =             L 1 0 0 0 0 0 L 2 0 0 0 0 0 0 H 1 0 0 0 H 3 0 H 2 0 0 0 H 4 0 0 H 1 0 -H 3 0 0 0 H 2 0 -H 4 0 0 0 0 L 1 0 0 0 0 0 L 2 0 0 0 0 -H 1,x 0 H 3,x 0 0 0 -H 2,x 0 H 4,x 0 0 H 1,x 0 0 0 H 3,x 0 H 2,x 0 0 0 H 4,x             , (5.10)
with L i and H i being respectively the Lagrange polynomials of the rst order and the Hermite polynomials of the third order. {s} is the vector of nodal displacements: {s} = {{s 1 } {s 2 }}, .11) Using the given denitions, eq. (5.8) can be expressed as a function of the nodal displacements:

s i = u t,i v c,i w c,i θ x,i θ y,i θ z,i . . ( 5
W int = {s} T [K el ] {s} , (5.12)
with [K el ] being the stiness matrix of the element. W ext is made by the contribution of inertial forces and external loads and it can be rewritten as a function of nodal displacements:

W ext = W inertial + W loads = -{s} T [M el ] {s} + {s} T {F el } , (5.13)
where [M el ] is the mass matrix of the element and {F el } is the vector of nodal loads. Equation (5.6) can be simplied using expressions (5.12) and (5.13). That allows to obtain the nal spatial discretization of the elasticity problem for a single element of the structure:

[M el ] {s} + [K el ] {s} = {F el } . (5.14)
The procedure to obtain the mass and stiness matrices and their elements will not be described here and the reader can refer to [START_REF] De Frias Lopez | A 3D nite beam element for the modelling of composite wind turbine wings[END_REF]. The linear system (5.14) can be extended to the whole structure by properly assembling the two matrices and the load vector, hence:

[M ] d + [K] {d} = {F } , (5.15)
where {d} is the vector of nodal displacements for the entire structure. Eq. (5.15) is discretized in the time domain using the Newmark method [START_REF] Newmark | A method of computation for structural dynamics[END_REF]:

     d(t + ∆t) = d(t) + δt (1 -γ) d(t) + γ d(t + ∆t) , d(t + ∆t) = d(t) + ∆t ḋ(t) + ∆t 2 2 (1 -2β) d(t) + 2β d(t + ∆t) , (5.16)
with 0 β 1 2 and 0 γ 1. If γ = 0.5 the method is second-order accurate regardless of β. If β > 0.25 the method is unconditionally stable, therefore no limits of δt are imposed by stability issues. For the present work γ = β = 0.5 is chosen.

Substituting equation (5.15) into (5.16), the fully discretized system of equation for the elasticity problem is obtained:

∆t 2 [K] [M ] [M ] + ∆t 2 2 [K] 0    {d k+1 } ḋk+1    =    [M ] ḋk -∆t 2 [K] {d k } + ∆t {F k } [M ] ({d k } + ∆t ḋk ) + ∆t 2 2 {F k }    , (5.17)
where

d k = d(t k ) and d k+1 = d(t k+1 ).
This model has been implemented and validated for the case of uniform loads. Indeed, a comparison of the simulated results has been done by colleagues at Inria with the analytic theory of beams [START_REF] Pezzano | Aeroelastic modelling of a wind turbine[END_REF].

1.2 Fluid-Structure Coupling in the Three-Dimensional Solver

Two methodologies exist to deal with numerical simulation of the uid-structure interaction [START_REF] Hou | Numerical methods for uid-structure interaction -A review[END_REF]: the monolithic and the partitioned approaches. The monolithic approach consists in a total coupling of the uid with the elastic body, leading to a single system of equations, whose solution describes the entire domain of the problem. On the other hand, the partitioned approach separates the elasticity problem from the Navier-Stokes equation and solves two distinct systems While the monolithic approach can achieve better accuracy and is well suited for immersed boundary method, the partitioned one is more exible and easier to implement. Indeed, it requires little modications to both the structure and the uid codes. Moreover, splitting the domains allows to use dierent meshes for the two, not only in term of computational cells, but also in term of topology, such as coupling a 3D nite volume method with a 1D nite element beam.

It was chosen to exploit the simplicity of the beam structure model and a partitioned scheme has been employed. That means, that, as represented in gure 5.2, during each time iteration, there is an exchange of information between the Navier-Stokes algorithm and the nite element model. The structural code needs indeed to know the aerodynamic forces to apply to the blade to evaluate the displacements. Those are then communicated back to the uid solver to update the geometry and the body velocity. Moreover, the computation of the level-set function needs to be updated given the modication of the geometry. This implies a signicant increase of the computational time, especially for complex geometries, as the blade, where a high number of Lagrangian markers is needed to properly describe it.

To summarize, a 6 degrees of freedom per node nite element beam model is employed with: linear axial deformation, bending, based on Euler-Bernoulli theory, linear torsion, based on Saint-Venant hypothesis, anisotropic stress-strain relationship for composite materials.

Indeed, a french company named Tensyl, specialist of composite structure and mechanical systems have been consulted for our wind blade model (APX48). They are able to compute highdelity structural simulations of wind blade deformations, which is not our aim in the present work. They explained that all the coupling phenomena are not signicant for this blade. Therefore, the terms AT , F T and LT , that represent respectively the axial-torsion, the ap-torsion and the lag-torsion coupling are not considered in our beam model. The mass distribution, axial and bending stinesses used in the structural model are those of the APX48 wind blade.

Rotating Framework for Wind Blade Modelling

A rst diculty was to convert the geometry of the blade from a standardized CAD format into an array of properly organized Lagrangian markers that will be the input for both the Octree Navier-Stokes solver and the nite element code presented above. We had available the design of the wind blade only in a format describing the blade with non-uniform rational b-splines.

A second obstacle has also been added with the presence of the hub of the rotor. The levelset function is so computed for both obstacles. The modications operated in the uid and structure models will now be described before presenting some preliminary results obtained in this conguration.

Boundary Conditions and Rotating Models

The rotating frame uid model will now be explained. Considering an inertial frame (X, Y, Z) and a rotating one (X , Y , Z ), the relationship between the velocity of a point expressed in the two references is:

V i = V r + ω ∧ r, (5.18)
with r the position in a rotating reference, ω is the rotational speed, subscript i means the inertial frame of reference and r means the rotating one. For the present work, the y-axis of the domain is taken as the rotation axis of the wind blade:

ω =     0 Ω 0     . (5.19)
It is therefore straightforward to obtain:

∇.V i = ∇.V r ∇ 2 .V i = ∇ 2 .V r (5.20)
In the case of a uid system, the acceleration is equal to the Lagrangian derivative:

a = DV Dt , (5.21)
and the relationship between the acceleration in the two references can be written:

a i = a r + 2ω ∧ V r + ω ∧ (ω ∧ r ) + dω dt ∧ r , (5.22)
where r' is the position in the rotating frame, ω ∧ (ω ∧ r ) represents the centrifugal acceleration and the term 2ω ∧ V r is the Coriolis acceleration. Since our blade is submitted to a constant revolution rate, and using the denition of eq. 5.21, equation 5.22 becomes:

DV i Dt = DV r Dt + ω ∧ (ω ∧ r ) + 2ω ∧ V r .
(5.23)

The Navier-Stokes momentum equation, in the Lagrangian form is:

DV i Dt = - 1 ρ ∇p + ν∇ 2 V i .
(5.24)

Substituting expressions 5.20 and 5.24 into 5.23 and considering the continuity equation, the incompressible Navier-Stokes equations for a rotating reference frame are obtained:

   ∇ • V r = 0 DV r Dt = - 1 ρ ∇p + ν∇ 2 V r -ω ∧ (ω ∧ r ) -2ω ∧ V r (5.25)
The two non inertial terms act as volume forces and they can be accounted in the prediction step of the Chorin algorithm. Using the expression of ω, the centrifugal and Coriolis accelerations are:

ω ∧ (ω ∧ r ) = ω ∧     Ωz' 0 -Ωx'     = -Ω 2     x' 0 z'     and (5.26) 2ω ∧ V r = 2Ω     w' 0 -u'     .
(

The boundary conditions used for the velocity eld have also been modied to simulate the rotating blade compared to previous simulations done for code validation. Indeed, the following conditions have been used for the velocity:

At the inlet (y-direction), Dirichlet boundary conditions with u y = U ∞ , u x = -Ωz and u z = Ωx are imposed.

At the outlet, homogeneous Neumann conditions are imposed.

In the lateral x and z directions, Dirichlet boundary conditions are used, exactly the same as for the inlet.

Considering now the structural model of the blade, the Coriolis and centrifugal forces need to be included in the nite element model described above. The method used is widely described in the Master thesis of Pezzano [START_REF] Pezzano | Aeroelastic modelling of a wind turbine[END_REF] who worked at MEMPHIS team. It was chosen to keep the structural model as simple as possible. Therefore, the rotation is assumed to be only with respect to the z-axis, which means that tilting and rotation movement of the hub are neglected. Moreover, the dynamic eects caused by other rotational deformation speeds are also neglected and gyroscopic torques are not considered. The work of the Coriolis force and of the centrifugal force are written thanks to their denitions. The structural equations are solved in order to obtain the matrix of the shape functions of the complete beam model. From equation 5.15, the space discretized equations of elasticity are then obtained by taking into account the newly added contributions:

[M ] d + [C] ḋ + [K] {d} = {F ext } + {F cf } = {F } , (5.28)
with [C], the Coriolis matrix and {F cf } is the contribution of the centrifugal force. Applying the Newmark method with γ = β = 0.5 allows to get the linear system for the solution of the dynamic structural problem:

∆t 2 [K] [M ] + ∆t 2 [C] [M ] + ∆t 2 2 [K] ∆t 2 2 [C]    {d k+1 } ḋk+1    =    [M ] -∆t 2 [C] ḋk -∆t 2 [K] {d k } + ∆t {F k } [M ] ({d k } + ∆t ḋk ) + ∆t 2 2 {F k }    .
(5.29)

Preliminary Results on a Rotating Wind Blade

Once the implementation of the rotating framework has been performed, all the features developed during the present work were ready to simulate a rotating blade with operating conditions corresponding to those of the real wind turbine. A thought has been done to dene the computational domain and the mesh to use in order to perform the simulations. We are indeed dealing with a blade whose length is 23.5 meters. A compromise has been done on the computational domain size in order not to reach a high number of degrees of freedom while having a correct resolution on the blade. Moreover, it is needed to be far enough from the boundary conditions. Otherwise, abnormal behaviours of the ow around the blade were observed. The boundary conditions have indeed a key impact on the numerical tool developed. The computational domain chosen is [-30, 30] x [-15, 45] x [-30, 30] and it is centred at the middle of the blade hub and the blade extends in the negative z-direction. The incoming wind arrives from the negative y-direction and y-axis corresponds to the rotation axis of the wind blade. This can be seen in gure 5.3.

The grid chosen is also a compromise to reduce the computational time. The minimum size of the grid is ∆x = 0.0293 and ve dierent levels of renements are employed. The objective being to simulate the uid-structure interaction occurring during the rotation of the blade, it was necessary to reduce the number of degrees of freedom. As explained briey previously, simulating with the structural model of the blade implies to compute the level-set function and its gradients at each time step which is very signicantly time consuming. Indeed, the number of Lagrangian markers necessary to properly describe the blade geometry is around 200 000. Moreover, the gradient of the total level-set (hub+blade) should be evaluated for the computation of the wall law, whereas for the aerodynamic forces computation, the use of the gradient of the blade level-set is necessary. Therefore, a compromise has been done on the minimum grid size and on the number of cells in order to use a reasonable number of computational hours. The total number of cells of the chosen grid is around 28 million and the octree mesh is displayed in gure 5.4.

All the features developed during the present work have been used for this preliminary simulation with the wall law, the turbulence model, the rotating framework. This section indicates preliminary results since, now, only a uniform wind speed is imposed at inow. So, this case does not constitute a real eld condition of the blade. However, a real operating condition has been simulated in the sense that an incoming wind of 8 m.s -1 with a rotational speed of 21.6 rpm has been used, which corresponds to one working condition of the wind turbine. The viscosity has been imposed as the one of the air with a usual temperature. As a way to validate the results, the power extracted by the simulated wind turbine is recorded. The gure 5.5 shows the streamwise velocity of the ow after 3.5 revolutions of the wind blade. The contour lines reveal the tip vortex phenomenon that is developing in the wake of the blade. When looking at the streamwise velocity in x, y and z-plane, no small structures can be seen and especially, in the z-plane, the wake behind the airfoil is not very extended. Indeed, as the Eulerian scheme to treat the advection term is only rst order accurate, signicant numerical dissipation occurs during the simulation. Moreover, the used mesh does not allow a very accurate simulation of the ow. The streamwise velocity in the y-plane shows the rotation movement of the blade. In the x-plane, the ow pattern of the tip vortices can be seen in the wake. Simulation of 22 s of ow (that corresponds to around 8 revolutions of the blade) has been performed on Occigen with 600 processors. The corresponding computational time was 338 hours. The power extracted by the wind turbine is computed during this simulation as a way to compare with the real operating data. The computation is done thanks to equation (1.4) that denes the power as a function of the torque of the three blades. Data from the wind turbine manufacturer inform us that the expected power extracted by this wind turbine with the previously detailed operating condition is 251 kW. After 8 revolutions of the blade, the simulated value begin to be stabilized and we obtained a power of around 114 kW as can be seen in gure 5.6. It can be concluded that a good order of magnitude is obtained. However, the simulation performed did not allow to get a good accuracy for this result. The possible improvements that could be done are the use of a second-order accuracy for the penalization and to rene in the area of the wind blade. Both possibilities would certainly imply a better accuracy in the computation of the aerodynamic forces. We are indeed simulating a very complex case. Also, even with the high numerical dissipation, the ow pattern with the tip-vortex phenomenon seems to represent correctly the ow around a wind blade.

Adaptive Mesh Renement Process

An adaptive mesh renement (AMR) process has been implemented in the code in order to optimize the number of cells while controlling the accuracy. The aim is indeed to rene in the most interesting areas. For our present work, the areas of highest interest in the ows studied are basically the wakes of the dierent objects tested. We could thus be able to rene the grid on the zone that follows the wake and coarsen elsewhere. An adequate criterion should be dened in order to comply with this objective. Several criteria have been tested for the case of the ow around a sphere. This ow has been chosen since it deals with three-dimensional eects and this case does not require much computational time to perform tests. The dierent criteria tested were: based on the q-criterion, based on the vorticity, based on the L 2 norm of velocity gradients normalized by the L 2 norm of velocity. For example, in 2D, the criterion is:

∂U x ∂x 2 + ∂U x ∂y 2 + ∂U y ∂x 2 + ∂U y ∂y 2 U 2 x + U 2 y .
The tests performed with a criterion for the AMR based on the q-criterion showed inadequate behaviours. Indeed, the q-criterion allows to highlight the vortices. Therefore, a grid adaptation based on it does not permit to follow the whole wake but only the vortices. The two other criteria show similar results and allow to adapt the grid along the wake as requested. The advantage of the criterion based on velocity gradients normalized by the norm of velocity is the fact that the threshold to indicate for coarsening and for renement is less dependant of the studied ow when using the normalization by the velocity. Indeed, depending on the case to simulate (dierent obstacles, dierent Reynolds numbers,...), the parameters for adaptation will vary. The chosen criterion shows to be very less dependant to these ow congurations and the values to indicate for adaptation can be dene with very few tests. AMR process with superimposition of the mesh at the top the two-dimensional solver with two dierent ows as a proof of concept. Indeed, for 2D cases, the use of grid adaptation is less pertinent as in 3D since the gain in term of cell numbers has less impact, moreover in light of the computational time needed for the AMR process.

Figure 5.7 shows the grid adaptation for the case of the ow around a circular cylinder at Re = 200. The same AMR process has been performed for the ow around a NACA0012 airfoil at Re = 1000. As an indication, the same values for the coarsening and renement parameters have been used for both ow congurations. The grid is actually rened in the wakes but the coarsening can be improved. It has been chosen to perform the grid adaptation only for the 3 ner level of trees in order not to have too close level jumps. No further work has been pursued to nd the right parameters for the two-dimensional case since we wanted to spend more eorts on the three-dimensional case. We validated the method on these two-dimensional cases before working in the 3D solver.

Three-Dimensional Cases

For the computation of the gradients in 3D case, the same method as used for the gradient and divergence operators in the Navier-Stokes equations, explained in the chapter 2, is employed. A particular attention has been brought to nd adequate values for the parameters to rene and coarsen the grid. It has also been chosen to perform the AMR for 3D cases only with the two ner levels of octree to avoid several level jumps in the same area that could imply a signicant loss of accuracy. The AMR has been tested rstly for the case of the ow around a sphere at Re = 5000. This quite high Reynolds number has been chosen in order to have a turbulent ow, more representative to the simulation of the rotating blade. The grid adaptation obtained for this case can be seen in gures 5.8 to 5.10 and we can observe that a ne optimization in the grid distribution is reached. Moreover, no loss of accuracy seems to occur at the level jumps.

The AMR process has then been tested for the case of the rotating blade which is our case of highest interest. A global quite coarse mesh has been used for the tests. It was noticed that the values to use for the grid adaptation needed to be changed compared to the case of the sphere. Indeed, the both ows are very dierent, especially since, for the blade, the rotating framework is employed. Figure 5.11 shows the AMR process for the ow around the rotating blade with a low Reynolds Number (around 1000). It can be observe that a good grid adaptation is reached. However, when increasing the Reynolds number and simulating with a real operating condition of the wind turbine, the mesh adaptation doesn't t properly with the wake of the blade. It shows that the criteria need to be slightly adapted. This work will be extended in order to perform future simulations of the rotating blade with AMR.

Rotating blade submitted to a gust

As a rst approach, a xed mesh has been used and a gust-like perturbation has been submitted to the rotating blade. The experimental gust events highlighted in the previous chapter show that the mean duration of an event is around 30 seconds. But we were not able with the remaining time (human and computational) to simulate such a long period. Therefore, we performed a simulation with a fast increase of the velocity to reproduce a gust-like behaviour. The preliminary simulation described above was used as an initial state, with a wind speed of 8 m.s -1 as boundary condition. The gust is then imposed as a linear increase of the wind speed. This linear increase lasts 0.5 s and the maximal wind speed of 13 m.s -1 is then reached. Gust events highlighted during the experimental work have been used to have ideas of reasonable values of the wind speed increase during a real gust. The boundary conditions in the lateral x and y directions have been adapted compared to those used for the preliminary simulation. As a reminder, Dirichlet boundary conditions were rstly imposed with u y = U ∞ , u x = -Ωz and u z = Ωx. But, as we have now an increase of U ∞ along time, it seemed more reasonable to impose Neumann boundary condition in the y direction for the 4 lateral boundaries of the computational domain.

The structural model of the blade has been included in the simulation as a rst approach, as well as the wall law and the turbulence model. The objective would be to determine the eects of the gust on the structural behaviour of the blade. An extension of the present work will indeed be performed in order to obtain a Reduced Order Model (ROM) of the structural model of the blade. Both the experimental data of pressure collected on the wind blade and the pressure data from this simulation will be used to build the Proper Orthogonal Decomposition (POD) modes of the structural model of the wind turbine. ROMs could allow to reduce significantly the computational time needed for the uid-structure interactions. The pressure data is so recorded at each Lagrangian marker that describes the blade geometry and at each time step. But the computational cost of the structural model with that mesh and ow conguration was way too important to be performed, even with a gust-like event faster than real gusts. As a comparison, with the same computational time, and in the case of a gust event, the simulation without structural model computes around 1 s whereas with the structural model only around 0.075 s was simulated. This was therefore not possible to use the structural model with that conguration. As previously explained, the cost of the uid-structure interaction is not only due to the structural model but also to the computation of the level-set function and its gradients.

The structural model has so been removed to perform the simulation. When the streamwise velocity were looked up during the gust-like event, an incorrect behaviour of the ow is observed. This can be seen in gure 5.12. At the beginning of the wind speed increase, the behaviour stays similar as previously showed, while having an increase of velocity in front of the blade. But after a few times, the ow does not behave correctly anymore. This could be due either to the change in the boundary conditions, compared to the case without gust, that might be unsuitable for this complex ow. Indeed, the solver developed is very sensitive to the choice of the boundary conditions. An other possibility could be the fact that the increase of incoming wind is too fast for the solver.

Finally, even if we are not ready to simulate a gust event, the preliminary simulations show promising results with an octree solver which is only rst-order accurate. The ow with the wall law and the turbulence model seems to reproduce the behaviour of the ow around a wind blade. The numerical tool has been validated in the previous chapter with several non-trivial test cases which is a promising result. The Adaptive Mesh renement process needs slight adjustments of the used criteria before being employed in the simulations, but it could allow to reduce the numerical diusion in the wake of the blade. We also have at our disposal a high number of experimental data with dierent kind of values as the wind speeds at dierent heights, that allows to have real boundary conditions, pressure data and deformations on the wind blade. This is of high interest to perform comparisons with several dierent simulations.

Conclusion

In this thesis we have proposed a numerical method to solve the incompressible Navier-Stokes equations on octree meshes. Two and three dimensional solvers have so been developed in that context with nite volume method. The application eld of this work is the study of the ow around wind turbines. Therefore, as dealing with high Reynolds numbers, Large Eddy Simulations with a subgrid turbulence model has been implemented. We were aware that octree grids don't allow to properly dene the anisotropic boundary layer around the wind blade. Indeed, an octree mesh can only rene by decreasing its size by two, which is not enough for the boundary layer. Consequently, the choice has been done to implement in the solver a wall law to model the smallest scales close to the blade. The structural behaviour of a blade plays a signicant part in the ow around it, so a uid-structure interaction has been developed with a beam model of the blade. The structural model, as well as the geometry of a real wind blade have been used.

The existing numerical models used for preliminary design of wind turbines are rst described. These tools, that are often partly based on the Blade Element Momentum theory have the advantage of being fast to compute. But, the geometry of the blade should be known with the polar curves of the airfoils, the chords, the twist angles... It is therefore not possible to perform simulations of a new blade or of a blade with aerodynamic appendices, without a preliminary work. On the other hand, Computational Fluid Dynamics (CFD) allow the computation of new geometries. The grid to use is a key parameter for CFD since it needs to be ne enough to accurately model the ow, while having a signicant importance on the computational time. Numerical schemes with a discretization based on octree grids have been investigated and these meshes show an interesting compromise between computational time and accuracy for our application to ow around wind turbines. It is indeed well adapted to Adaptive Mesh Renement process, that allows to rene in the areas of high interest in the computational domain while coarsen elsewhere.

In the present work, the Navier-Stokes equations are solved with a penalization method which is well adapted to Cartesian or octree grids. In this framework of penalized Navier-Stokes equations, a two-dimensional solver based on a octree discretization is developed. Finite volume method is mainly used in the thesis. The numerical scheme developed to treat the Laplace operator is based on Discrete Duality Finite Volume methods with a Diamond dual mesh introduced to compute the discrete gradient at the face centre. For the advection term, a semi-Lagrangian scheme has been developed that have several advantages as it avoids the computation of the convective term and is unconditionally stable. All the operators developed are validated with the use of analytical functions and their consistency is proved. These methods are then extended to the three-dimensional case and the operators are again validated with the use of analytical functions that prove their consistency. As previously explained, with an application of the solver being the modelling of the ow around wind blades, wall functions for high Reynolds number ow have been implemented together with Large Eddy Simulations. A Vreman subgrid model has been used and implemented in the context of octree grids.

As the models proposed in the present work are original and have the objective of being a compromise between accuracy and computational time, a wide validation was essential. The consistency and accuracy of the whole solver in 2D has been computed thanks to the analytical solution of the Taylor Green vortex. A second order accuracy is reached in space and a rst order in time. Data from literature has then been used to perform comparisons with the solver developed. The ow around a circular cylinder has been studied at several Reynolds numbers. The drag coecient and the Strouhal number have been compared to the literature and our results show to be in line with the literature. The ow around a NACA0012 airfoil has also been studied to perform simulations with a more complex geometry. The lift coecient computed with our solver has thus been validated by comparing with literature. While performing simulations with the three-dimensional solver, a stability issue occurred with some grids. It has been realized that it was due to the semi-Lagrangian scheme. As we were not able to x that issue, a rst-order Eulerian scheme has been implemented instead to solve the advection terms of Navier-Stokes equations. The 3D solver has then been widely tested. The ow around a sphere has been simulated at low Reynolds number and the drag coecient has been validated by comparing with data from literature. The ow around a cylinder at Reynolds 3900 has then been widely studied with a comparison of wake proles simulated with the developed solver and with experimental data. After this ow conguration at a medium Reynolds number, a simulation with the ow around the cylinder at higher Reynolds, 140 000, has been performed. This ow simulation was performed with LES and a viscous wall function. This case has been compared with experimental data from literature and also with a well-known Cartesian code. The grid used for both simulations is compared and arguments are given concerning the computational time. The results obtained are in-line with literature, which allowed to validated the numerical tool developed.

In the framework of a H2020 European Project (AeroGust), we had the opportunity to perform an experimental work in order to collect in-service wind turbine data. The aim of the project is to study the eects of gust on a wind turbine blade. Specications for this study were to get measurements on a wind blade where we know geometrical and structural data. Wind measurements are also needed in order to represent the wind as it arrives at the wind turbine. Valorem uses usually meteorological masts to perform wind assessment studies for wind farm development as prescribes in the IEC 61400-12-1 standard [1]. As the wind during a gust is characterized by high variations in wind velocities in a few seconds as explained in the IEC standard, the acquisition frequency of wind data should be small enough. In order to enable estimation of the vertical wind prole (strong gradients can exist in the vicinity of the ground), anemometers at 5 dierent heights have been installed. Concerning the blade data, Valorem has not the expertise to install sensors on a strong experimental case which is an operating wind blade. So, a complete technological solution was looked up with regards to sensors supply, installation, maintenance and communication system. A global requirement was to minimize the impact of measurement tools on the blade as the wind turbine should be working well during and after the experimental phase. In order to have experimental data showing the aero-elastic behaviour of a wind blade, the measurement of blade deformations will allow to observe its structural behaviour. To observe the aerodynamic load on the wind blade, the measurement of pressure of air will be of signicant interest. The technological solution of choice uses optical bre with Fiber Bragg Gratings (FBG) sensors. They measure the wavelength change induced by deformation or other data. The bres have been stuck along the blade surface and connected to an inquiry system located in the hub of the wind turbine. The same sampling of 1 s is done.

Filtering and analysing of the data has been performed but it is still ongoing to identify the inuence of wind speed on the deformations measured. Several gust events have been identied and we have now at our disposal several test cases with many measurements of wind speed, deformations and pressure at several locations.

The structural model of the wind blade that has been instrumented, has been implemented in the developed solver with a simple beam model. Indeed, one application of the present work is the study of the eects of gusts on blade deformations. A weak uid-structure coupling is performed in the numerical tool, with two dierent solvers, one for the uid and one for the structure, that exchange only the necessary informations. In order to reduce the size of the computational domain, the choice has been done to perform simulations of one rotating blade with an inertial domain. Adaptations have thus been done in the uid and structural models to take into account the Coriolis and centrifugal forces in the computations. Simulations have so been performed in that framework with a uniform incoming wind speed and a real operating condition of the wind blade. A lot of numerical dissipation occurred in that simulation, which is due to the rst-order accuracy of the solver and to the fact that the grid used was not very ne in the wake of the blade. An Adaptive Mesh Renement process has been implemented in the solver and has been investigated with two and three dimensional cases. Promising results have been obtained, but a slight adaptation of the criterion should be done to nd the appropriate thresholds of the criterion to rene and coarsen the grid for the rotating blade at high Reynolds ow. A simulation has then been performed with a gust-like event whose behaviour has been dened from the experimental data. The computational cost of the uid-structure coupling is very high, therefore, a perspective of the present work is the use of Reduced Order Models with Proper Orthogonal Decomposition to model the structural beam model of the wind blade.

Perspectives

As we saw that the uid-structure coupling is very costly, the use of Reduced Order Models (ROMs) would be of high interest. We could use both the experimental data and data from CFD to compute Proper Orthogonal Decomposition (POD) modes. The pressure at every Lagrangian marker that denes the geometry of the blade is indeed recorded during the simulation which allows to have at our disposal a high number of CFD data. Thanks to the experimental work, we also have signicant experimental values with a lot of wind conditions. In a more general way, modelling of the future performance of a wind plant (taking into account the terrain, the wake of other wind turbines) is of high interest for Valorem. Since the computational time is a key parameter for the engineers that have to deal with a high number of wind resource assessments, ROMs could be an interesting eld of research, coupled with the present numerical tool and with the experimental data collected.

The solver developed in the thesis constitutes a basic architecture with several features (LES, wall law, AMR, uid-structure interaction). This base can now be improved quite easily, especially by having a second-order accurate Eulerian scheme to compute the advection term of the Navier-Stokes equations. A method is currently investigated in the Memphis team of INRIA. The velocity eld is known at the cell-centre of each cell. The gradients of velocity are then computed on each cell centre with a Moving Least Square interpolation (that method is described in chapter 3, section 4). The velocity uxes are thus rebuilt at every face centre thanks to the gradients and velocity eld from cell centres. This is done in an upwind framework. Another upgrade would consist in getting second-order accuracy in time. For example, an Adams-Bashforth scheme could be implemented for the convective term and a Cranck-Nicholson method could be used to compute the viscous term. That would yield second-order accuracy in time for the prediction step. These upgrades have been implemented in the solver by colleagues in the Memphis team in a few weeks. New simulations have been performed and allowed to observe the small structures of the turbulence in the wake of cylinder at Reynolds 140 000 with the same octree mesh as presented in chapter 3. This can be seen in the gure below. That shows that the solver developed during the present thesis constitutes an interesting reliable numerical tool and continuous improvements can be conducted.

List of Figures

Résumé détaillé

Résolution des équations de Navier-Stokes sur maillage octree : vers une application à la modélisation d'une pale d'éolienne Contexte de la recherche Pendant la Conférence de Paris de 2015 sur les changements climatiques, presque tous les états du monde ont trouvé un accord sur la réduction des changements climatiques. Un des objectifs annoncé est de limiter le réchauement climatique à moins de 2°C comparé à l'ère pré-industrielle. L'accord vise également à atteindre une émission de gaz à eet de serre due à l'activité humaine à un niveau nul. Certes, la conférence ne planie pas d'actions concrètes pour les pays qui ont validé l'accord, et il n'est pas prévu de sanction dans l'hypothèse où un état n'eectuerait pas d'eort susant. Tout de même, cet accord montre que le monde est en train de réaliser l'importance pour tous de limiter les changements climatiques. An d'atteindre les objectifs achés, le développement des énergies renouvelables va jouer un rôle fondamental et il semble raisonnable de penser qu'un scénario gagnant serait le développement d'un mix énergétique avec du solaire, de l'éolien, de l'hydraulique, de la méthanisation ... À l'intérieur de ce mix, le secteur éolien présente plusieurs avantages qui le distingue des autres moyens de production. La ressource en vent est en eet locale et son coût baisse de plus en plus grâce aux améliorations technologiques. De plus, contrairement à l'énergie solaire, le vent est disponible à toute heure. La France a aussi l'avantage de disposer de nombreux couloirs de vent distincts et cela permet de disposer presque toujours au moins sur une partie du territoire de la ressource éolienne. Le développement de la lière éolienne oshore pourrait également permettre d'accélérer en France la production éolienne. L'éolien oshore est en eet en train de bénécier d'avancées technologiques qui permettent de baisser les coûts d'installation petit à petit. Actuellement, le coût au MWh de l'éolien oshore est plus bas que le coût produit par l'EPR de Flamanville (France) en cours de construction.

Ce développement de la lière éolienne entraîne ainsi de nombreux avancements technologiques. La gure 5.14 illustre ce phénomène avec l'évolution de la taille des éoliennes au cours du temps. Grâce aux longues pales, les éoliennes sont en eet capables de produire même lorsque la ressource en vent est limitée. En ce qui concerne l'entreprise Valorem, elle intervient tout le long de la chaîne de valeur des énergies renouvelables, du développement des parcs jusqu'à l'exploitation et la maintenance en passant par la construction. Le développement d'un parc éolien consiste à identier des zones d'implantations potentielles avec les contraintes associées et obtenir les diérentes autorisations administratives nécessaires. Lors de cette étape, plusieurs études sont eectuées an d'estimer la ressource en vent disponible et la future production d'énergie (en fonction des modèles d'éoliennes et de leur disposition sur le terrain). Des études environnementales sont faites pour estimer les impacts du parc. Également, des études de raccordement électriques sont intégrées lors de cette phase de développement pour préparer le raccordement du parc au réseau. Valorem possède plusieurs parcs (éoliens majoritairement,

Démarche adoptée

Dans un premier temps, l'aérodynamique des éolienne est étudiée an de comprendre les paramètres importants à calculer. Les outils numériques simpliés permettant d'eectuer des calculs rapidement pour estimer la puissance extraite par une éolienne sont étudiés, en particulier la méthode Blade Element Momentum, qui est beaucoup utilisée par les industriels du domaine. Les modèles numériques de CFD (Computational Fluid Dynamics) sont ensuite considérés. Comme il existe de très nombreuses méthodes, nous nous sommes intéressées à celles dont l'application porte sur la modélisation des éoliennes, et également aux schémas numériques qui utilisent une discrétisation sur des maillages octree. Ces types de maillages ont en eet été étudiées et montrent des avantages intéressant pour notre application. Lors de cette thèse, des méthodes numériques ont été développées an de résoudre les équations de Navier-Stokes dans un contexte de frontières immergées. Le choix de cette méthode est justié dans le chapitre 2. Un code de calcul utilisant des méthodes type volumes nis avec une discrétisation sur des maillages octree a été développé dans cette thèse. Dans un premier temps un code en deux dimensions a été développé. Les méthodes d'interpolations avec les stencils utilisés sont décrits. Chaque opérateur qui permet de résoudre les équations de Navier-Stokes est décrit et validé dans un premier temps, en utilisant des fonctions analytiques. Ces schémas numériques sur maillage octree sont ensuite modiés an de résoudre le problème en trois dimensions. Comme notre domaine d'application est l'écoulement de l'air autour des pales d'éoliennes, des grands nombres de Reynolds doivent être simulés. Une méthode Large Eddy Simulation (LES) avec un modèle de turbulence pour les petites structures de type Vreman a été implémentée. An de modéliser les petites échelles proche de la paroi, une méthode dite des wall functions a été développée. En eet, l'utilisation de maillages type octree ne permet pas de modéliser correctement la couche limite proche paroi.

Une validation détaillée des code 2D et 3D a ensuite été eectuée. Cette étape est une partie essentielle du travail de la thèse. Le code de calcul 2D a été testé pour le cas de l'écoulement autour d'un cylindre 2D pour plusieurs régimes d'écoulement et ensuite pour le cas de l'écoulement autour d'un prol aérodynamique. Grâce aux données provenant de la littérature, les résultats des simulations ont pu être comparés avec diérentes géométries et diérents nombre de Reynolds. Le code de calcul 3D a ensuite également été largement contrôlé en procédant par étapes. Dans un premier temps, le cas de l'écoulement autour d'une sphère à faible nombre de Reynolds (500) a été considéré. Ensuite, des simulations ont été eectuées pour simuler l'écoulement autour d'un cylindre à Re = 3900. Le modèle de turbulence, ainsi qu'une loi de paroi visqueuse ont été intégrés aux simulations. Enn, un haut nombre de Reynolds de 140 000 a été simulé. Les données de la littérature sur ce cas ont également été utilisées pour discuter sur les diérences entre un code de calcul d'ordre élevé sur maillage Cartésien et la présente thèse.

Un travail expérimental a été conduit pendant la thèse. En eet, très peu de données expérimentales sont disponibles avec à la fois la connaissance de la géométrie de la pale, le modèle de structure et des mesures sur la pale. Et aucun travail disponible ne permet d'avoir accès, dans le même temps, à la donnée de la vitesse du vent qui arrive sur la pale. Le contexte de la participation de Valorem au projet européen AeroGust 1 a permis de collecter des données expérimentales sur une éolienne en fonctionnement. Une réexion a été menée an de déterminer le cahier des charges permettant de mener à bien cette étude expérimentale. Plusieurs solutions techniques ont été investiguées an d'être capables de mesurer le vent et également pour obtenir les données sur la pale et le choix nal est présenté. L'installation sur le terrain a été organisée et eectuée. Les données ont ensuite été récupérés pendant plusieurs mois et analysées.

Une application du présent travail est la modélisation de l'écoulement autour des éoliennes. Comme nous avons à disposition des données expérimentales de déformations d'une pale d'éolienne en fonctionnement, le modèle de structure de cette pale a été implémenté dans le code octree. Un couplage uide-structure est développé. Comme l'écoulement autour d'une éolienne complète aurait nécessité un coût de calcul trop important par rapport aux ressources disponibles, une seule pale a été étudiée dans nos simulations. De plus, un système de référence tournant a été mis en place pour focaliser le domaine sur la pale. Une fonctionnalité intéressante a également été développée avec un processus d'adaptation automatique du maillage (AMR). Cela permet d'optimiser le maillage en étant le plus précis possible dans les zones d'intérêt (comme le sillage) tout en diminuant le nombre de mailles loin de la pale. Les réelles conditions de fonctionnement de l'éolienne ont ainsi été simulées.

Principaux résultats obtenus

Les équations de Navier-Stokes sont écrites dans le contexte d'une frontière immergée avec une méthode de pénalisation de l'obstacle. On considère un domaine Ω qui comprend à la fois un uide et un (ou plusieurs) obstacle(s). La frontière extérieure de ce domaine est notée Γ. Cependant, durant la phase de validation du solver pour le cas 3D, un problème de stabilité s'est produit pour certains maillages octree. L'origine de ce problème de stabilité s'est avéré venir du schéma semi-Lagrangien. Il a donc été décidé de changer de méthode et de développer un schéma Eulérien upwind à l'ordre 1 à la place.

Une validation principale a été eectuée avec le cas de l'écoulement autour d'un cylindre à Reynolds 140 000. Le modèle de turbulence ainsi que la loi de paroi ont été intégrées dans les simulations et une comparaison est faite entre le code de calcul développée dans la thèse et un code de calcul Cartésien, d'ordre élevé qui a été développé au préalable de cette thèse. Ainsi, le code octree permet d'avoir un nombre de degrés de libertés signicativement plus faible que dans le cas Cartésien. Le temps de calcul n'est toutefois pas seulement dépendant du nombre de mailles, le conditionnement des matrices a un impact également et elles sont bien mieux conditionnées dans le cas d'un maillage Cartésien. Une première visualisation sur la gure 5.15 permet d'observer que les petites échelles de la turbulence sont captées par le code Cartésien mais pas le code octree.

Des données expérimentales provenant de la littérature ont également été utilisées an de comparer les prols de vitesse dans le sillage du cylindre avec les deux codes de calcul. Une première simulation a été eectuée avec les deux codes an de passer la phase transitoire de démarrage des calculs. Les résultats visibles sur la gure 5.16 ont ensuite été obtenus avec une moyenne sur 5 lachés tourbillonnaires. Cela permet de valider le code de calcul développé dans cette thèse.

Le procédé d'adaptation automatique du maillage a été développé pour plusieurs applications. Un eort supplémentaire est encore nécessaire pour adapter le critère utilisé pour le rendre performant pour les hauts nombres de Reynolds dans le système de référence tournant. Des simulations préliminaires ont été eectuées sur une pale en conditions réelles de fonctionnement. Malgré une forte dissipation numérique due à l'ordre peu élevé des schémas numériques, la forme de l'écoulement observé est en adéquation avec la réalité. Le code de calcul développé constitue ainsi une base de travail avec de nombreuses fonctionnalités. Les schémas numériques développés pourront par la suite faire l'objet d'améliorations an d'augmenter leur précision. La campagne de mesures expérimentales nous permet de disposer d'une importante base de données qui pourra être utilisée pour réaliser des futures simulations de conditions réelles de fonctionnement de l'éolienne et permettre de calibrer les modèles numériques.

Figure 1 :

 1 Figure 1: Evolution of the size of wind turbines over time. Source www.windeurope.org

Figure 1 . 2 :

 12 Figure 1.2: C Pλ curve for the Euros EU120 wind blade. Source www.euros.de.

Figure 1 . 3 :

 13 Figure 1.3: Example of power curve for some usual wind turbine.

Figure 2 . 2 :

 22 Figure 2.2: Examples of quadtree grids used for validation

Figure 2 . 3 :

 23 Figure 2.3: Example of mesh conguration in 2D with stencils for Laplacian operator

Figure 2 . 4 :

 24 Figure 2.4: Example of error distribution between a Laplace operator discretization and an analytical expression and numerical solution with a quadtree grid conguration corresponding to tree level 7

Figure 2 . 5 :

 25 Figure 2.5: Example of mesh conguration in 2D with stencils for semi-Lagrangian scheme

Figure 2 . 6 :

 26 Figure 2.6: Example of octree grid conguration used for validation of the solver

Figure 2 . 8 :

 28 Figure 2.8: Turbulent boundary layer proles for dierent Reynolds numbers 4.2 Near Wall Modelling Despite the growth in computer performances, wall resolved Large Eddy Simulation still remains unfeasible to model the ow around wind turbines. It has been estimated by Choi et al.[START_REF] Choi | Grid-point requirements for large eddy simulation: Chapman's estimates revisited[END_REF] that the number of grid points required by a wall resolved LES is of the order of Re

Figure 2 . 9 :

 29 Figure 2.9: Reichardt's wall law. Found in [2].

  Figure 3.1: Velocity eld of the Taylor-Green vortex on a quadtree grid (tree level = 7) after 0.1 s of simulation

Figure 3 . 2 :

 32 Figure 3.2: Quadtree mesh around a two-dimensional cylinder (minimum tree level = 7)

Figure 3 . 3 :Figure 3 . 4 :Figure 3 . 5 :

 333435 Figure 3.3: Drag coecient history at Re = 550 with the numerical solution of reference and the computed results with 3 dierent grid sizes

Figure 3 . 6 :

 36 Figure 3.6: Isobars and contour lines of z-component of vorticity obtained at t = 100 and Re = 200

Figure 3 . 7 :

 37 Figure 3.7: Drag and Lift coecients history at Re = 200

  (a) Case with semi-Lagrangian prediction (b) Case with Eulerian prediction

Figure 3 . 10 :Figure 3 . 11 :

 310311 Figure 3.10: Simulation of a blade in a rotating frame, y-Velocity at plane z = -5

3 Figure 3 . 12 :

 3312 Figure 3.12: Time history of drag coecients for the ow at Re = 500 with 3 dierent grid renements and 2 dierent CFL conditions

5

  level jumps, domain Ω = [-8; 16] × [-12; 12] × [-12; 12] respectively in the x, y and z directions, minimum size of the grid ∆x = 1.172.10 -2 , number of points on the diameter of the cylinder = 85, total number of cells = 138 011 392.

Figure 3 . 14 :

 314 Figure 3.14: Pressure coecient around the cylinder at Re = 3900. Angle of π corresponds to the wake of the cylinder, whereas angle of 0 is the incoming ow.

Figure 3 . 15 :

 315 Figure 3.15: Wake prole of streamwise velocity at dierent positions obtained by an average over 9 vortex sheddings after a preliminary simulation for the ow past a cylinder at Re = 3900.

- 8 ,

 8 16] × [-8, 8] × [0, 4] . The boundary conditions are: Dirichlet boundary condition for the velocity at the inlet (x-direction) homogeneous Neumann condition for the velocity at outlet periodic conditions in the lateral z-directions that correspond to the cylinder length homogeneous Neumann condition for the velocity in the lateral y-direction.

Figure 3 . 16 :

 316 Figure 3.16: Comparison of x-Velocity eld obtained with the Cartesian and the octree codes for the ow around a cylinder at Re = 140000

  (a) Streamwise velocity at x = 1 (b) Normal velocity at x = 1 (c) Streamwise velocity at x = 3

Figure 3 . 17 :Figure 3 . 18 :

 317318 Figure 3.17: Wake prole for the ow past a cylinder at Re = 140000 at dierent locations obtained by an average over 5 vortex shedding after a preliminary simulation

Figure 4 . 1 :

 41 Figure 4.1: Computational domain with d = rotor diameter + a few meters

Figure 4 . 2 :

 42 Figure 4.2: Topographic view of the wind plant obtained with Google Earth ® and wind rose corre- sponding to the site.

Figure 4 .

 4 4 shows an illustration of the instrumentation of the blade with this technology.

Figure 4 . 3 :

 43 Figure 4.3: Fiber Bragg Grating principle. n corresponds to the index of refraction and Λ is the period of n variation. Found on www.fbgs.com.

Figure 4 . 4 :

 44 Figure 4.4: Optical bre technology

Figure 4 . 6 :

 46 Figure 4.6: Position of pressure sensors on the blade

Figure 4 .

 4 Figure 4.8 shows a photograph of the whole mast after its installation. Figure4.9 presents the sensors of the met mast with anemometers and wind vanes.

Figure 4 .

 4 Figure 4.8 shows a photograph of the whole mast after its installation. Figure4.9 presents the sensors of the met mast with anemometers and wind vanes.

Figure 4 . 7 :

 47 Figure 4.7: Final position of the met mast and wind rose of the area.

Figure 4 . 8 :

 48 Figure 4.8: Photo of the met mast after its installation

Figure 4 . 9 :

 49 Figure 4.9: Photo of the sensors on the met mast

Figure 4 . 10 :

 410 Figure 4.10: Aerial work platform positioning

Figure 4 . 11 :

 411 Figure 4.11: Photo of the pressure side of the wind blade after installation of sensors

  (a) Weak wind data (b) Gust event

Figure 4 . 13 :

 413 Figure 4.13: Time history of deformation data from sensor located on the leading edge at 7m of the blade length, of wind speed data reconstructed at the wind turbine and frequency (with a normalized scale) during a gust and with weak wind.

  (a) Deformations at pressure side with weak wind (b) Deformations at pressure side during gusts (c) Deformations at trailing edge with weak wind (d) Deformations at trailing edge during gusts

Figure 4 . 14 :

 414 Figure 4.14: Deformations on the blade at several location as a function of blade position with weak wind (a and c) or during gust events (b and d)

Figure 5 . 2 :

 52 Figure 5.2: Partitioned scheme implemented in the 3D solver

Figure 5 . 3 :

 53 Figure 5.3: Sketch of the computational domain for the ow past a rotating blade.

Figure 5 . 4 :

 54 Figure 5.4: Display of the octree mesh used for the simulation of the ow past a rotating blade.

Figure 5 . 5 :

 55 Figure 5.5: Display of the ow around a rotating blade with U ∞ = 8 m.s -1 and Ω = 21.6 rpm after 3.5 revolutions of the blade.

Figure 5 . 7 :

 57 Figure 5.7: Velocity magnitude for the ow past a circular cylinder in 2D at Re = 200 showing the

Figure 5 . 8 :Figure 5 . 9 :Figure 5 . 10 :

 5859510 Figure 5.8: Beginning of the simulation with AMR process for the case of the ow around a sphere at Re = 5000

Figure 5 . 11 :Figure 5 . 12 :

 511512 Figure 5.11: AMR process for the rotating blade with an operating condition at U ∞ = 6 m.s -1 , Ω = 16.2 rpm and Re = 1000

Figure 5 . 13 :

 513 Figure 5.13: Streamwise velocity past a cylinder for the ow at Re=140 000.

1

  Evolution of the size of wind turbines over time. Source www.windeurope.org . . 1.1 Sketch of the aerodynamic forces exerted by the ow around a wind blade airfoil 1.2 C Pλ curve for the Euros EU120 wind blade. Source www.euros.de. . . . . . . . 1.3 Example of power curve for some usual wind turbine. . . . . . . . . . . . . . . . . 1.4 Axial stream tube around a wind turbine [41] . . . . . . . . . . . . . . . . . . . . 1.5 Blade element theory [41]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.6 Comparison of predicted power for BEM and GAD models with experimental data of NREL UAE phase 6 (experimental process described in [28]). . . . . . . . 1.7 Example of quadtree discretization and corresponding tree representation. Found in [74]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.8 Sketch of the computational domain . . . . . . . . . . . . . . . . . . . . . . . . . 2.1 Example of mesh conguration in 2D with stencils . . . . . . . . . . . . . . . . . 2.2 Examples of quadtree grids used for validation . . . . . . . . . . . . . . . . . . . . 2.3 Example of mesh conguration in 2D with stencils for Laplacian operator . . . . 2.4 Example of error distribution between a Laplace operator discretization and an analytical expression and numerical solution with a quadtree grid conguration corresponding to tree level 7 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Example of mesh conguration in 2D with stencils for semi-Lagrangian scheme . 2.6 Example of octree grid conguration used for validation of the solver . . . . . . . 2.7 Example of mesh conguration in 3D with stencils for Laplace operator . . . . . 2.8 Turbulent boundary layer proles for dierent Reynolds numbers . . . . . . . . . 2.9 Reichardt's wall law. Found in [2]. . . . . . . . . . . . . . . . . . . . . . . . . . . 2.10 Implementation of the wall law. Found in [26]. . . . . . . . . . . . . . . . . . . .3.1 Velocity eld of the Taylor-Green vortex on a quadtree grid (tree level = 7) after 0.1 s of simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Quadtree mesh around a two-dimensional cylinder (minimum tree level = 7) . . . 3.3 Drag coecient history at Re = 550 with the numerical solution of reference and the computed results with 3 dierent grid sizes . . . . . . . . . . . . . . . . . . . 3.4 Isobars and contour lines of z-component of vorticity obtained at t = 100 and Re = 100 for the ow around a circular cylinder . . . . . . . . . . . . . . . . . . . 3.5 Drag and Lift coecients history at Re = 100 . . . . . . . . . . . . . . . . . . . . 3.6 Isobars and contour lines of z-component of vorticity obtained at t = 100 and Re = 200 for the ow around a circular cylinder . . . . . . . . . . . . . . . . . . . 3.7 Drag and Lift coecients history at Re = 200 . . . . . . . . . . . . . . . . . . . . 3.8 Isobars and contour lines of z-component of vorticity obtained at t = 100 and Re = 600 for the ow around a circular cylinder . . . . . . . . . . . . . . . . . . . 3.9 Drag and Lift coecients history at Re = 600 . . . . . . . . . . . . . . . . . . . .

Figure 5 . 14 :

 514 Figure 5.14: Évolution de la taille des éolienne au cours du temps. Source www.windeurope.org

  -u) dans Ω, ∇.u = 0 dans Ω, u(x, 0) = u 0 (x) dans Ω, u = u f (x, t) sur Γ f , où p est la pression, u représente le champs de vitesse. La variable ρ est la masse volumique du uide, ν sa viscosité cinématique et χ B est la fonction caractéristique :χ B = 0 in theuid domain 1 in the solid domain Ainsi, un code de calcul a été développé pour résoudre ces équations dans le contexte parallèle d'un maillage de type octree. La méthode des volumes nis a été principalement utilisée dans cette thèse. Le schéma numérique développé pour traiter l'opérateur de Laplace est basé sur une méthode DDFV (Discrete Duality Finite Volume) où un maillage dual en forme de diamant a été introduit et qui permet de calculer les gradients discrets au centre des faces de chaque cellule. Dans le cas du terme d'advection, un schéma semi-Lagrangien a été implémenté dans un premier temps. Cette méthode a en eet l'avantage d'éviter le calcul du terme de convection et de ne pas être soumis à la contrainte de stabilité où à la condition CFL (Courant-Friedrich-Lévy).

Figure 5 . 15 :

 515 Figure 5.15: Comparaison de la vitesse selon l'axe x obtenue avec les codes Cartesien et octree pour le cas de l'écoulement autour d'un cylindre à Re = 140000

  Incompressible Navier-Stokes equations are solved with a nite volume method. Pressure and velocity variables are collocated. A block structured mesh is used with general curvilinear coordinates. Large Eddy Simulation (LES) is performed by applying a low-pass lter to the Navier-Stokes equations, resulting in ltered velocity eld. A sub-grid model is then employed to model the scales below the grid size. This solver is well documented and is open-source under the MIT License.The solver Overow CFD has been studied since a lot of features are available in this research code developed by NASA Ames Research Center. Chaderjian described in

Table 2 .

 2 1: Norm of the error and order of the divergence operator for 2D Cartesian grid

Table 2 .

 2

2: Norm of the error and order of the divergence operator for 2D quadtree grid

Table 2 .

 2 3: Norm of the error and order of the Laplace operator for 2D Cartesian grid

	Tree level	L ∞	Order	L 2		Order
	5	2.837.10 -5		1.989.10 -5
	6	7.574.10 -6	1.91	5.394.10 -6	1.88
	7	1.973.10 -6	1.94	1.410.10 -6	1.94
	8	5.095.10 -7	1.95	3.564.10 -7	1.98
	9	1.252.10 -7	1.98	9.033.10 -8	1.98
	10	3.150.10 -8	1.99	2.272.10 -8	1.99
	Tree level -Number of cells	L ∞	Order		L 2	Order
	5 -1792	1.028.10 -3		4.254.10 -4
	6 -7168	2.522.10 -4	2.03	1.085.10 -4	1.97
	7 -28672	6.272.10 -5	1.93	2.750.10 -5	1.98
	8 -114688	1.566.10 -5	1.97	6.928.10 -6	1.99

9 -458752 3.917.10 -6 1.97 1.739.10 -6 1.99 10 -1835008 9.658.10 -7 2.01 4.306.10 -7 2.01 Table 2.4: Norm of the error and order of the Laplace operator for 2D quadtree grid

Table 2 .

 2 5: Norm of the error and order of the semi-Lagrangian for 2D quadtree grid

	Tree level -Number of cells	L ∞	Order	L 2	Order
	5 -1792	1.459.10 -3		3.875.10 -4
	6 -7168	7.025.10 -4	1.05	9.856.10 -5	1.98
	7 -28672	3.453.10 -4	1.02	2.486.10 -5	1.99
	8 -114688	1.707.10 -4	1.02	6.241.10 -6	1.99
	9 -458752	8.488.10 -5	1.01	1.564.10 -6	2.00
	10 -1835008	4.232.10 -5	1.00	3.914.10 -7	2.00

Table 2 .

 2 7: Norm of the error and order of the Laplace operator for 3D octree grid

	Tree level -Number of cells	L ∞	Order	L 2	Order
	5 -11264	5.764.10 -3		5.186.10 -3
	6 -90112	3.415.10 -4	2.03	1.339.10 -4	1.95
	7 -720896	8.492.10 -5	2.01	3.409.10 -5	1.97
	8 -5767168	2.121.10 -5	2.00	8.608.10 -5	1.99
	9 -46137344	5.312.10 -6	2.00	2.164.10 -6	1.99
	Tree level -Number of cells	L ∞	Order	L 2	Order
	5 -11264	2.681.10 -3		1.599.10 -3
	6 -90112	6.686.10 -4	2.00	4.702.10 -4	1.77
	7 -720896	1.671.10 -4	2.00	1.255.10 -4	1.91
	8 -5767168	4.178.10 -5	2.00	3.232.10 -5	1.96
	9 -46137344	5.312.10 -6	2.00	2.164.10 -6	1.98

table 3

 3

	.2.

Table 3 .

 3 2: Norm of the spatial error and orders for 3D octree grid

Table 3 .

 3 3: Two-dimensional grid characteristics and computational time for a simulation of 6 s

	tree level	Minimum size	Number of cells on the diameter	Global number of cells	Number of processors	Computational time
	7	1.465.10 -2 68	442 516	72	36 min
	8	7.324.10 -3 136	1 764 136	120	6 hours 56 min
	9	3.662.10 -3 273	7 091 050	480	35 hours 49 min

  Flow at Re = 200

	Re	Authors	S t	C D
	100	Braza et al. [12]	0.1600 1.3600
		Williamson (experimental) [91] 0.1640	-
		Henderson [39]	0.1664 1.3500
		He et al. [37]	0.1670 1.3528
		Bergmann [8]	0.1660 1.4100
		Present work	0.1736 1.4193

Table 3 .

 3 4: Comparison of Strouhal number and drag coecient for the ow past a circular cylinder at Re = 100

Table 3 .

 3 5: Comparison of Strouhal number and drag coecient for the ow past a circular cylinder at Re = 200can indeed be seen with the Strouhal number, whose value is higher as conrmed in table 3.6. The results obtained with the developed quadtree solver overestimated slightly the data from literature in this ow regime. The drag and lift coecient history of this simulation can be seen in gure 3.9. This gure was used to compute the averaged drag coecient and the Strouhal number of this ow.

	Re	Authors	S t	C D
	600	Henderson [39]	0.2294 1.4682
		He et al. [37]	0.2306 1.4641
		Henderson [38] (bi-dimensional tendency curve)	-	1.47
		Present work	0.2351 1.5096

Table 3 .

 3

6: Comparison of Strouhal number and drag coecient for the ow past a circular cylinder at Re = 600

  .7 where the averaged lift coecient at dierent angles of attack is reported. The quadtree solver shows accurate results in this ow conguration.

	Re	Authors	α = 10°α = 15°1
	000	Mittal et al [58]	0.42	-
		Khalid et al [43]	0.39	0.58
		Suzuki et al [86]	-	0.72
		Kurtulus [45]	0.42	0.70
		Present work	0.38	0.66

Table 3 .

 3 7: Comparison of lift coecient for the ow past a NACA0012 airfoil at Re = 1000 with dierent angles of attack

Table 3 .

 3 8: Three-dimensional grid characteristics and computational time for a simulation of 100 s

	5	4.688.10 -2 21	590 892	24	8 hours
	6	2.344.10 -2 42	4 711 596	144	25 hours
	7	1.172.10 -2 85	37 597 820	960	36 hours

Table 3 .

 3 9: Comparison of averaged drag coecient for the ow past a sphere at Re = 500

	C D

Table 3 .

 3 10: Comparison of averaged drag coecient for the ow past a cylinder at Re = 3900

	C D

Table 3 .

 3

	28

11: Comparison of averaged drag coecient (over 9 vortex sheddings) for the ow past a cylinder at Re = 140000

Table 4 .

 4 gusts f weak wind Amp gusts Amp weak wind 2: Mean values of frequency (noted f ) and amplitude (noted Amp) of deformations during gusts or with weak wind for each sensor

	Pressure side 2m	0.325	0.161	23.2	43.1
	Pressure side 7m	0.332	0.163	53.3	84.0
	Pressure side 14m	0.326	0.163	48.0	78.1
	Suction side 2m	0.321	0.162	23.4	44.6
	Suction side 7m	0.337	0.183	52.8	70.4
	Suction side 14m	0.323	0.163	48.1	72.9
	Trailing edge 2m	0.270	0.172	31.1	74.5
	Trailing edge 7m	0.266	0.171	87.2	180.1
	Leading edge 2m	0.267	0.160	28.2	30.0
	Leading edge 7m	0.264	0.160	86.0	146.4

www.aerogust.eu

Optimad, bitpit https://github.com/optimad/bitpit

Extension of the Methods to the Three-Dimensional Case

www.cines.fr

http://www.idris.fr/

Remerciements -Aknowledgments

Remerciements for the French National Centre for Scientic Research (CNRS).

Ministry of Higher Education and Research) and acquired by GENCI (Grand Equipement National de Calcul Intensif) allowed to perform major simulations of the present work. Finally, the supercomputer Turing, operated by IDRIS (Institute for Development and Resources in Intensive Scientic Computing) has been used for the simulations with

Collection and Analysis of Operating Wind Turbine Data

Wind shear and gusts are both wind situations that create the largest loads on wind turbines. Indeed, the magnitude of a gust can be suciently high to increase the loads up to the point of fracture. Catastrophic failures are reported in the literature for both wind turbine blades and towers (see for example [START_REF] Chou | Failure analysis and risk management of a collapsed large wind turbine tower[END_REF] and [START_REF] Chou | Failure analysis of wind turbine blade under critical wind loads[END_REF]).

Thus, the task of collection of experimental data on an in-service wind turbine has been planned in the context of the European funded project AEROGUST (More informations on www.aerogust.eu). Interests will be rst to have real data and to use it to better understand the eects of wind and more precisely of gusts, on wind blades. A second interest is to use experimental data to calibrate our numerical schemes in the high-delity CFD code described above. Indeed, very little experimental data of wind turbines is accessible which includes the geometry, structural model and measurements on blade (pressure, strain). Further, none of the available data correlate these measurements with upstream gusts coming on to the turbines, to the best of our knowledge. As an extension for future application, wind turbine aeroelastic ROMs (Reduced Order Models) could be created from a combination of high delity CFD and data collected from an in-service turbine. The initial plan (before the beginning of AEROGUST) was to obtain experimental data from wind tunnel measurements but the costs were prohibitive. As Valorem Group owns and operates wind farms and has many contacts with French wind farm owners, the idea of obtaining experimental data from an in-service wind turbine was deemed more viable.

The motivations of this experimental work are thus to:

Analyse eects of wind on blade loads.

Compare experimental data with those obtained with our CFD code.

A requirement for the choice of the wind turbine to measure, is that the geometry and the structural model of the blades are known. This is needed to meaningfully compare experimental measurements and numerical simulations. As Valorem only has those data for one wind blade model (APX48 model), which is not very common, the choice of the wind turbine to measure was limited. The Valorem Group, on the other hand, doesn't own this kind of wind turbine but has already worked in partnership with a wind plant owner having Jeumont J48 wind turbines (a model with APX48 blades). APX48 was therefore the blade model of choice for this project.

On the following sections, the work done by Valorem and INRIA to specify the data to be The adaptive mesh renement process has been implemented in the 2D and 3D solvers. This stage is performed after the prediction step so that the zero divergence condition is guaranteed at the end of the grid adaptation. Indeed, once the grid has been adapted, the predicted velocity eld computed at the prediction step should be updated for the new grid conguration. This is done in the present work in a simple way. Thanks to the bitpit library, dierent functions exist so that each cell of the new grid conguration has information toward the old conguration. Three dierent possibilities exist for each cell:

if the cell has the same size as in the former conguration, the value contained in the former cell is applied to the "new" one;

if the cell has been rened since the former state, the value contained in the former bigger cell is applied;

if the cell has been coarsened, a basic mean is performed with all the values contained in the smaller cells constituting this new bigger cell. For example, in 2D, the coarsened cell will have the value of 1 4 of the sum of the 4 former smaller cells.

The level-set function and its gradient should also be computed again after each grid adaptation, which consist in the most time-consuming step of the AMR process. The dierent results obtained with several test-cases will now be described.

Two-Dimensional Cases

The computation of the velocity gradients is thus needed to perform the grid adaptation. The same method as used for the divergence and gradient operators in the Navier-Stokes solver as described in chapter 2 for two-dimensional case is used. The AMR process has been tested in Résumé :

Le sujet de la thèse est le développement d'un outil numérique qui permet de modéliser l'écoulement autour des pales d'éoliennes. Nous nous sommes intéressés à la résolution des équations de Navier-Stokes en incompressible sur des maillages de type octree où les échelles plus petites en proche parois ont été modélisées par la méthode dite des wall functions. Abstract :

The subject of the thesis is the development of a numerical tool that allows to model the ow around wind blades. We are interested in the solving of incompressible Navier-Stokes equations on octree grids, where the smallest scales close to the wall have been modelled by the use of the so-called Wall Functions. An automatic Adaptive Mesh Renement (AMR) process has been developed in order to rene the mesh in the areas where the vorticity is higher. The structural model of a real wind blade has also been implemented and coupled with the uid model. Indeed, an application of the numerical tool is the study of the eects of wind gusts on blades. An experimental work has been conducted with an in-service wind turbine with the measurement of wind speed upstream. This data will allow to calibrate and validate the numerical models developed in the thesis.

Key words :

Incompressible ows; Octree grid; Immersed boundary; Fluid-Structure Interaction; Wind Turbine