Keywords: Fluid-structure interactions, SPH-FE coupling, Aquaplaning------

Avant toute chose je tiens à remercier la DGA d'avoir refusée de financer le premier sujet de thèse à laquelle j'avais postulé. Sans ce refus, je ne serais sans doute pas en train d'écrire ces mots. De plus j'aurais passé les trois dernières à travailler sur une méthode maillée, ce qui me semble aujourd'hui une hérésie!! Ce constat étant fait, je tiens à remercier en premier une personnalité atypique, ce cher Bibenbum! Merci à toi de m'avoir accompagner durant ces 3 années, et de m'avoir insufflé une meilleur façon d'avancer! Connaissant mon amour de l'écriture et de la langue anglaise, tu as su trouver les mots pour me fournir le courage de finir ce manuscrit dans les temps.

Après Bibenbum, je tiens à remercier mon directeur de thèse David Le Touzé et mon premier responsable entreprise Christian Mignot pour avoir m'avoir offert l'opportunité d'entreprendre cette thèse. Ce qui m'a permis de jongler entre les intérêts d'un laboratoire publique et d'une entreprise industrielle. Je remercie également David pour sa bienveillance et son expertise durant ces 3 années. Tes remarques et conseils ont toujours été pertinentes. Merci pour ton implication pour cette thèse qui s'éloigne beaucoup du milieu industriel avec lequel travaille principalement ton équipe. Je remercie Christian et Bibendum de m'avoir fournis les moyens de travailler confortablement, sans avoir d'autres préoccupation que de me concentrer à la science sans recherche de financement additionnels. Toujours chez Michelin, je souhaite remercier Yohan Le Chenadec et Violaine Todoroff pour leurs regards critiques et leurs conseils avisés pour garder en permanence en mémoire la visée finale de ces travaux. Parmi mon encadrement, je souhaite remercier plus particulièrement mon encadrant de thèse Guillaume Oger pour son aide tout au long de ces années au LHEEA. Au delà de tes conseils, de ta patience et de ton expérience c'est surtout de ta bonne humeur quotidienne et de ton optimisme dont je te suis reconnaissant. Sans cela, il aurait bien plus difficile de finir ce travail de recherche. Tes remarques constructives et tes coups de boosts m'ont clairement aidé à arriver jusqu'au bout. Je conseille d'ailleurs aux futurs doctorants de toujours se retrouver à côté de Guillaume lors d'une conférence, tant pour les remarques physiques qu'humoristiques.

Après mon encadrement, je remercie également Antonio Souto-Iglesias et Laurent Gornet d'avoir suivi mon travail durant ces trois années mais également d'avoir accepté de faire partie de mon jury de thèse. Je vous remercie d'avoir pris le temps et le courage de rentrer dans le sujet. Je souhaite également remercier Thomas Rung, Andrea Colagrossi et Christophe Berthon d'avoir accepter de participer à ce jury. J'ai une pensée particulière pour les jurés français, venus travailler un jour férié.

Je souhaite ensuite remercier l'ensemble du laboratoire, enseignants, permanents, postdoctorants, doctorants pour la bonne ambiance et pour cette première expérience dans le monde la recherche. Je n'oublie pas Anne, Sonia et Elodie pour leurs coups de mains administratifs! J'ai une pensée particulière pour mes anciens enseignants de l'option Océan, Lionel, Guillaume, Félicien et David qui m'ont donné envie de débuter dans la recherche. Je remercie également mes compagnons de recherche, les doctorants. Je souhaite bon courage à mon successeur Julien pour le chemin qui se dresse devant lui, même si je ne me fais pas de soucis avec l'encadrement dont il bénéficiera au laboratoire! J'ai une pensée particulière pour Laurent mon compagnon de recherche et de bureau. Je te remercie d'avoir partagé tes gâteaux et de m'avoir supporté durant ces années malgré mes chansons, l'odeur du sac de rugby et autres bizarreries. Merci à toi pour ces longues conversations scientifiques à refaire la SPH, ou encore pour ces belles rigolades en regardant les images de ce bon vieux Paul. J'ai également une pensée particulièrement aux Marie's, brune et blonde, du bureau d'à côté. Nous avons décidé tous les trois de se lancer dans cette aventure à la suite de l'option Océan sans vraiment savoir ce qui nous attendait vraiment. Merci à toi Marie, pour ces bonnes années de colocation à la rue Racine où il faisait bon vivre malgré les cinq étages à monter avec les courses. En dehors de tous ces camarades, je souhaite remercier mon entourage. Merci aux copains du rugby qui m'ont permis de bien me défouler sur le prés et de bien se vider la tête après le travail. Merci aux rugbeux de Centrale et de l'ARCCN pour tous les bons moments. Merci aux copains Kervégan pour toute l'intelligence et la délicatesse dont vous savez faire preuve. Je sais que je peux toujours compter sur vous. Merci également aux copains du VSN avec qui j'ai connu joie et peine pendant ces différentes saisons de rugby, finissant à chaque fois en finale! Une pensée particulièrement pour deux jeunes gens en or, le jeune Mabille et le fougueux Magique Labache. Ne changez pas! Je souhaite également remercier ma famille pour leurs soutiens et leurs encouragements. Enfin je remercie grandement Louloute qui m'a accompagnée pendant tout ce trajet. Je te remercie également pour ta patience et ton écoute attentive pour mes problèmes de travail! Merci pour ta folie et ta générosité au quotidien.

Pour finir, je remercie l'ensemble des personnes venues assister ma soutenance de thèse. Elles ont eu le courage de se lever tôt un jour férié pour observer l'aboutissement de ce travail de recherche. J'ai une pensée particulièrement à Malo, Rachel, Quentin, Dédé, Gigi et Mathilde qui sont venus sachant pertinemment qu'ils ne comprendraient pas grand chose. Je pense également à Niska, notre jeune compagnon à quatre pattes qui est très certainement le premier chien à avoir assister à une soutenance au laboratoire. Merci à toi d'avoir été sage et de ne pas avoir mangé la moquette! CONTENTS LIST OF FIGURES 5.8 Time history of E V (dashed lines) and E C (lines) of the water for different spatial resolutions in the fluid, Newmark scheme with α s = -0.3. . . . . . . . 5.9 Time history of E f luid/body (lines) and E body/f luid (dashed lines) for different spatial resolutions in the fluid, Newmark scheme with α s = -0.3. . . . . . . . 5.10 Time history of ∆E F +S (lines) and E Outside (dashed lines) for different spatial resolutions in the fluid, Newmark scheme with α s = -0.3. . . . . . . . . . . . 5.11

List of Symbols

Motivations

This research fits into a partnership between Michelin and the LHEEA laboratory of Ecole Centrale de Nantes. Michelin is a leader maker in the highly competitive tire market and is therefore constantly seeking to improve its tire performances for security purposes but also for societal issues. For example, maintaining the grip performances for a worn tire allows to extend the tires life, reducing their environmental impact. The risk of the aquaplaning phenomenon is well known of all drivers. This phenomenon represents a challenging problem both for tire and car industries. It strongly influences the tire grip properties on wet roads, especially for road braking and holding. It corresponds actually to a loss of contact surface between the tire and the road. A part of the water is evacuated into the grooves while the other one is discharged in front of the tire, imposing a strong pressure on the tire. When a critical velocity is reached, vertical fluid forces become too large to keep the contact between the tire and the road, affecting the road holding. More precisely, the shape of the tire's contact patch, its grooves and sipes, push the water forwards and drain away part of the water builds up in front of the tire. The water still creeping in between the contact patch and the road surface is then channelled into the grooves where it is stored. Angled and transversal tire grooves are first of all designed to drain away to the side as much of the water as possible. The tire can thus break through the residual water film and restore direct contact with the road surface. These three water dispersal stages (Fig. 1.1) occur between the leading and the lagging edge of the contact patch and correspond to three different transition zones called (ALBERT, 1968):

• Hydrodynamic zone • Viscodynamic zone • Damp zone, intermittent residual film of water
In the hydrodynamic zone, the dispersal and drainage effects are predominant. The impact of the tire on the water at the front of the contact patch causes a rise in the water pressure (called hydrodynamic pressure). If this pressure becomes greater than the mean pressure of the tire on the road surface, the tire can no longer repel the water and it lifts off the road surface. The hydrodynamic pressure increases in proportion to the square of the speed, [START_REF] Horne | Phenomena of pneumatic tire hydroplaning[END_REF]. Since the 60's, a large range of researches have been conducted to estimate the aquaplaning speed with the influencing factors such as the inflation pressure, the wheel load, the tire footprint ratio, etc. (see [START_REF] Horne | Phenomena of pneumatic tire hydroplaning[END_REF], [START_REF] Gengenbach | Experimental investigation of tires on wet pavements[END_REF][START_REF] Horne | Recent studies to investigate effects of tire footprint aspect ratio on dynamic hydroplaning speed[END_REF]). A rounded shape of contact patch understandably ploughs more easily through the water than a rectangular shape. At any given depth of water, the narrower the tire, the smaller the flow of water to be channelled through the grooves. On the one hand, a narrow contact patch therefore reduces the volume of water to be channelled through the contact patch. On the other hand, since wide tires have to deal with a much greater flow of water, a much greater water dispersal capacity must be built into the design. Note that the aquaplaning mechanisms are clearly different between car and truck tires. Today, the aquaplaning speed of a high-performance car fitted with correctly-inflated new tires exceeds 100 km.h -1 .

The viscodynamic zone corresponds to the water storage in the tire grooves. The tire design at the front of the contact patch, in the hydrodynamic zone, reduces to about few microns to 0.1 mm the depth of water that creeps under the contact patch. The water remaining must be channelled into the tire grooves where it is stored. The rubber blocks in the contact patch do not roll: they are laid down on the road surface at the leading edge of the contact patch, somewhat comparable to the way we put our feet on the ground as we walk. When the blocks then leave the contact patch, they simply lift off. It could be said that tire blocks compress water almost vertically. It is this compression that drives the water towards the grooves. There is one condition to be complied with for tire-road contact to be restored: the water compressed must be able to escape to the edge of each block before it leaves the contact patch. If the water cannot escape quickly enough, the rubber block does not touch the road surface. It those cases a question arise about the fluid compressibility. Note that it seems impossible to measure experimentally such potential compressible effects. The higher the pressure on the film of water and the shorter the distance to the edges of the rubber blocks, the shorter the transfer time to the storage 1.1. MOTIVATIONS zone will be. This would appear to mean that the rubber blocks should be small and the grooves large. How does the viscosity influences the phenomenon? It seems reasonable to considered a turbulent flow around the tire, but what kind of flow do we have inside the grooves or upstream the contact patch? How does it evolves with the water height on the ground? A study of the time for water transfer to the storage zones has been done by [START_REF] Wiess | Influence of Pattern Void on Hydroplaning and Related Target Conflicts[END_REF]. Nevertheless, this study does not consider the eventual capture of air bubbles inside the grooves. Ventilation effects could strongly modify this time transfer.

In the damp zone, the previous dispersal of water is not enough to make the road completely dry underneath the tire. Traces of water remain on the road surface, just as traces of water remain on the surface of a receptacle that has just been emptied. The surface tension of water causes microscopic drops of water to remain on any surface that has been covered by water. This presence of water modifies the grip mechanisms such as the friction laws through some lubrication effects at microscopic scales. Molecular adhesion can no longer operate if contact between the rubber and the road surface is not perfectly clean and dry. Indentation is inoperative if the depth of water covers the tiny bumps in the road surface. In wet ground, maintaining grip therefore involves dispersion the water to restore dry contact between the tire and the road. The tire must therefore find a way of restoring dry contact, in spite of this film of residual water a few microns thick.

Even if we only consider the fluid problem, many other physical questions arise. For an example the strong de-pressurization, which is occurring inside the grooves, may lead a cavitation phenomenon. A better understanding of the aquaplaning generating mechanisms is therefore needed to design safer tires. The numerical tools seem extremely attractive to answer this issue, in view of predicting tire performances and complementing the (expensive) experimental tests. However, numerical models and experimental tests are complementary tools. Simulation may be able to provide additional informations which are still unavailable through experiments. On the other hand, experiments are still critical for validation purposes.

The aquaplaning phenomenon such as many violent FSI problems involve some nonnegligible coupling effects due to the presence of strong reciprocal interactions. The growth of computational power has enabled their modelling, resulting in the development of various numerical methods to model these coupled phenomena. In the last decades, mesh-based methods were preferred on this topic. The Finite Element (FE) method is classically used for structure modelling [START_REF] Bathe | Finite Element procedures[END_REF], while several methods can be used for fluid modelling such as Finite Volume (FV), Finite Difference (FD) or FE methods. The aquaplaning problem has been the topic of simulation works (see [START_REF] Cho | Numerical investigation of hydroplaning characteristics of three-dimensional patterned tire[END_REF][START_REF] Vincent | Augmented lagrangian and penalty methods for the simulation of two phase flows interacting with moving solids. application to hydroplaning flows interacting with real tire tread patterns[END_REF][START_REF] Kumar | Study of hydroplaning risk on rolling and sliding passenger car[END_REF]) emphasizing its complexity: fluid-structure interactions (FSI), highly deformable body, structure modelling with complex materials involved, contact with asphalt and the complexity of the resulting fluid flow (extremely complex interface, road dry up, ventilation, possible development of turbulence and cavitation). All these aspects lead to challenging problems for numerical simulations of both fluid and solid domains. Most of the literature related to the aquaplaning problem is focused on the use of mesh-based methods using between FD-FE [START_REF] Oh | Hydroplaning simulation for a straight-grooved tire by using FDM, FEM and an asymptotic method[END_REF], [START_REF] Kim | Hydroplaning simulation for tires using FEM, FDM and an asymptotic method[END_REF]), FE-FE (KOISHI et al. (2001)) and also FV-FE coupling strategies (OKANO andKOISHI (2001), NAKAJIMA et al. (2000)). Each of these methods possess its own specificities concerning the viscous model, the compressible fluid, etc. The presented aquaplaning con-CHAPTER 1. INTRODUCTION figurations usually considered high water height values, from 4 to 10 mm. It seems that these models have some difficulties to consider smaller heights. Note that during a strong storm, the typical maximum water height value is around 1 mm. It also appears that these mesh-based approaches have some difficulties to deal with the different complex interfaces of the phenomenon: the free-surface and the fluid-structure interfaces. This thesis aims at designing a coupling strategy able to handle complex free-surface fluid flows in the presence of deformable structures.

State of the art of coupling techniques

Two distinct coupling strategies can be cited: monolithic and partitioned approaches, [START_REF] Farhat | Robust and provably second-order explicit-explicit and implicit-explicit staggered time integrators for highly non-linear compressible fluidstructure interaction problems[END_REF]). In both formulations, the main idea is to correctly represent the coupling interface using admissible boundary conditions in both sub-domains. This means to define a coupling strategy with information transfers suited to the considered problem. The continuity of fluid and solid quantities (velocity and stress vector) must be ensured through the interface.

Monolithic coupling strategy

In the monolithic approach, fluid and solid governing equations are solved simultaneously or synchronously [START_REF] Michler | A monolithic approach to fluid-structure interaction[END_REF]. They represent a strong formulation for FSI resolutions. The continuity of variables through the coupling interface is naturally fulfilled. The main advantage of this method is to maintain (to a certain extent) the accuracy and the stability of the scheme. Such strategies are usually performed using a single software.

An example of monolithic coupling strategy is the one adopted by [START_REF] Deuff | Extrapolation au réel des mesures de pressions obtenues sur des cuves modèle réduit[END_REF] or [START_REF] Antoci | Numérical simulation of fluid structure interaction by SPH[END_REF]. They proposed a SPH-SPH coupling strategy to simulate complex FSI problems. The solid and fluid sub-domains are modelled using a SPH formulation. No interface detection is needed. [START_REF] Hwang | Development of a fully lagrangian MPS-based coupled method for simulation of fluid-structure interaction problems[END_REF] and [START_REF] Shao | A three-dimensional coupling method for fluidstructure interaction problems by using explicit MPS method and Hamiltonian MPS method[END_REF] developed the same strategy using the Moving Particle Semi-implicit (MPS) method, proposing a MPS-MPS coupling strategy for FSI problems. Interactions between solid and fluid particles are just computed using different properties for the two media. Nevertheless, the use of SPH or MPS methods for solid mechanics is limited to a very small range of applications, since these methods need some improvements to correctly model elastic bodies. Another technique is to use the Particle Finite Element Method (PFEM) from each sub-domain as proposed by [START_REF] Idelsohn | Interaction between an elastic structure and free-surface flows: experimental versus numerical comparisons using the PFEM[END_REF].

A monolithic approach is particularly difficult to achieve in most cases. For complex problems, its use is often prohibitive. Fluid and solid sub-domains have their own specificities, time scales, physics or computational considerations, etc. Indeed, implicit solution leads to preconditioning difficulties. On the contrary, explicit solution leads to prohibitive CPU times due to the presence of very small time steps in the structure. The use of different numerical methods is therefore often preferred, adapted for each medium. Different software are therefore needed to represent complex coupled phenomena (see DEGROOTE et al. ( 2010)).

STATE OF THE ART OF COUPLING TECHNIQUES

Partitioned coupling strategy

Partitioned couplings solve solid and fluid equations separately or alternately in time, FE-LIPPA et al. (1999). This is the most popular strategy used for FSI problems.

The boundary conditions are imposed asynchronously at the interface. The strategy adopted to perform information transfers is critical [START_REF] Farhat | Two efficient staggered algorithms for the serial and parallel solution of three dimensional nonlinear transient aeroelastic problems[END_REF]. It can be responsible for a shift in the temporal integration. Sometimes this shift is responsible for an increase in energy, inducing a modification of the convergence and accuracy properties. The stability of the coupling strategy is therefore not ensured. It depends on the parameters in each sub-domain, but also on the coupling algorithm itself. The coupling procedure needs the information transfers to be performed properly to maintain a sufficient accuracy in each sub-domain. Nevertheless, this weak coupling has the advantage of providing some degrees of freedom for both model components. Optimisations can be then performed independently on each solver. However, note that the same solver can be used for each sub-domain as done by [START_REF] Peseux | Hydrodynamic impact: Numerical and experimental investigations[END_REF], where a FE-FE strategy is used to represent deformable cone impacts.

In practice, the solid sub-domain is usually modelled using the Lagrangian description as we focus on the displacements or the deformations. Furthermore, the mesh-based FE method is generally preferred to particle method (such as SPH (GINGOLD and MON-AGHAN, 1977)) for its accuracy and stability properties. Conversely, different approaches are possible for the fluid sub-domain using Lagrangian, Eulerian or ALE descriptions, either mesh-based or particle methods. The choice of the fluid solver (and hence of the coupling strategy) depends on the flow physics to be modelled. Each of them presents advantages and drawbacks. The most popular schemes are presented here.

ALE -Lagrangian formulations

Contrary to the Eulerian description, the Lagrangian description can naturally account for the motion of the boundaries of the domain. Nevertheless, a Lagrangian description of the fluid would suffer from mesh distortions. The presence of large deformations of the fluid domain leads to the need of re-meshing procedures for the adaptation between fluid and solid meshes, which tend to increase the algorithmic complexity together with the computational costs. Using an Arbitrary Lagrangian Eulerian (ALE) description, fluid equations are not solved at fixed points but rather at some arbitrarily moving nodes, so that the fluid mesh can follow the solid deformations. Such tools have been successfully applied to many FSI problems [START_REF] Donea | An arbitrary Lagrangian-Eulerian Finite Element method for transient dynamic fluid-structure interactions[END_REF][START_REF] Donea | Arbitrary Lagrangian-Eulerian Methods[END_REF] or [START_REF] Ribet | Numerical modeling of the impact on water of a flexible structure by explicit Finite Element Method -comparisons with radioss numerical results and experiments[END_REF]. This strategy allows one to reduce the use of re-meshing tools but only partially. In this context, the Lagrangian FE method is usually preferred for the solid sub-domain, whereas FV or FE methods are often considered for the fluid.

Euler -Lagrangian formulations

Another possibility is to consider fixed fluid meshes crossed by a solid mesh during its motion and deformations, such as the Immersed Boundary Method (IBM) [START_REF] Peskin | The immersed boundary method[END_REF] or the Fictitious Domain Method [START_REF] Van Loon | Comparison of various fluid structure interaction methods for deformable bodies[END_REF]. The Eulerian description is usually adopted in computational fluid dynamics. The key point of such methods is the CHAPTER 1. INTRODUCTION link between Eulerian and Lagrangian variables. This type of procedure is very popular and gives good results for large solid deformations in confined flows.

Nevertheless, for the free surface flows in interaction with a structure tracking the moving free surface is a relatively difficult task. The free surface is usually described using interface methods such as Volume of Fluid (VOF, HIRT and NICHOLS (1981)) or Level-Set (OSHER and SETHIAN (1988)). These methods are difficult to handle (interface smearing) for complex interface problems, especially in presence of complex moving geometries as in the targeted aquaplaning applications.

Particle method -Lagrangian formulations

All the previously described methods used mesh-based descriptions with their related difficulties. Alternatively, particle methods can be adopted to avoid the complications linked to the mesh treatments, such as the MPS method [START_REF] Koishi | Hydroplaning simulation using fluid-structure interaction in LS-DYNA. LS-DYNA conference[END_REF]), the Finite Volume Particle Method (FVPM, QUINLAN et al. (2014)), or the Particle Finite Element Method (PFEM, IDELSOHN et al. (2008)). Obviously the SPH [START_REF] Gingold | Smoothed Particle Hydrodynamics: theory and application to non-spherical stars[END_REF]) method can also be employed. Each of these methods has its own specificities, but in all of them large deformations of the fluid sub-domain are therefore authorized. The boundary treatment remains a difficult problem due to the loss of consistency from the lack of neighbouring particles in the kernel support and the absence of direct discretization of the boundaries. However, meshless particle methods reduce the complexity of fluid-structure interface handling, especially in the presence of complex free surface flows and solid geometries.

The SPH method has been successfully applied to a large range of applications such as rigid body impacts [START_REF] Oger | Twodimensional SPH simulations of wedge water entries[END_REF], sloshing [START_REF] Souto-Iglesias | Liquid moment amplitude assessment in sloshing type problems with smooth particle hydrodynamics[END_REF], etc. The MPS method is similar but suffered from some accuracy and stability issues, [START_REF] Lee | Fluid-shell structure interaction analysis by coupled particle and finite element method[END_REF]. Nevertheless, it can provided correct results with several corrections, [START_REF] Khayyer | Enhancement of stability and accuracy of the moving particle semi-implicit method[END_REF]. On the contrary, [START_REF] Idelsohn | Interaction between an elastic structure and free-surface flows: experimental versus numerical comparisons using the PFEM[END_REF] proved the capacity of a coupling strategy between PFEM and FE methods to correctly model FSI problems. However, the whole domain needs to be discretized, leading to high computational costs. This is a problem for the targeted aquaplaning application where modelling the air phase is not necessarily needed.

Various SPH-FE coupling methods have already been proposed since the 90's (see [START_REF] Johnson | Linking of lagrangian particle methods to standard finite element methods for high velocity impact computations[END_REF][START_REF] Attaway | Coupling of Smooth Particle Hydrodynamics with the Finite Element method[END_REF]). Most popular strategies are based on a master-slave coupling, in which contact forces are computed to prevent the penetration of SPH particles into the FE mesh. Each method proposed has its own specificity. For instance, DE VUYST (2005) considered FE nodes as SPH particles. GROENENBOOM and CARTWRIGHT (2010) used a master-slave scheme to compute a penalty force to prevent interface penetration and to satisfy the momentum conservation. Furthermore, SPH-FE coupling can be used for FSI or solid impact problems [START_REF] Zhang | Coupling of Smoothed Particle Hydrodynamics and Finite Element method for impact dynamics simulation[END_REF][START_REF] Caleyron | Modeling of reinforced concrete trough SPH-FE coupling and its application to the simulation of a projectile's impact onto a slab[END_REF]. Nevertheless, no contact algorithm dedicated to avoid material interpenetration is needed when the SPH solver is accurate. [START_REF] Fourey | Développement d'une méthode de couplage fluide structure SPH Eléments Finis en vue de son application á l'hydrodynamique navale[END_REF][START_REF] Yang | SPH Simulation of Fluid Structure Interaction Problems with Application to Hovercraft[END_REF] showed that no special treatment is needed for the deformable body interface. They thus showed the interest of SPH-FE coupling for complex FSI problem modelling without contact algorithms for free surface flows using elastic structures. Both validated their coupling strategy 1.3. THESIS ORGANIZATION on several FSI benchmarks with increasing complexity, showing the capability of such implementation enabling to easily couple SPH with different FE software. FOUREY (2012) also analysed the SPH-FE coupling robustness according to various coupling parameters, highlighting the need for a sufficient dissipation on the high frequency loading to maintain stable calculations. More recently LI (2013), [START_REF] Li | A non-intrusive partitioned approach to couple smoothed particle hydrodynamics and finite element methods for transient fluid-structure interaction problems with large interface motion[END_REF] proposed a SPH-FE coupling strategy, then improved by [START_REF] Nunez Ramirez | A partitioned approach for the coupling of SPH and FE methods for transient nonlinear FSI problems with incompatible time-steps[END_REF], ensuring energy conservation through the fluid-structure interface. Nevertheless this technique is strongly intricated with the FE model, and therefore implies hard-coding within the solid solver.

Discussion

The present work focuses on a violent fluid-structure interaction problem involving complex free surface flows with deformable structures. These aspects lead to some difficulties when using traditional mesh-based methods, in which the fluid mesh follows the deformations of the structure. Large deformations of the fluid domain lead to the need for re-meshing tools and expensive computational costs. In this context, the interest of SPH for modelling efficiently the aquaplaning flow has been demonstrated by BARCAROLO et al. (2014a) using a rigid body assumption. The meshless and Lagrangian natures of SPH naturally avoid the problem of fluid/solid grid compatibility: the ground/fluid and the fluid/tire interfaces. It is perfectly suitable for such simulations involving large deformations in both fluid and solid sub-domains. Then Michelin, Ecole Centrale Nantes and NextFlow Software have recently tested the ability for the SPH solver to classify tires based on their surface structure geometries, without considering the gas phase. The other significant advantage of the SPH method, in this context, consists in its ability to be coupled with conventional Finite Element (FE) solvers, see [START_REF] Fourey | An efficient FSI coupling strategy between smoothed particle hydrodynamics and finite element methods[END_REF].

A partitioned coupling approach has been retained with the formalism introduced by FOUREY (2012). The solid part is modelled with a classical FE method. On the fluid side, the SPH method is chosen for its capability to handle the complex interfaces. The aim of this work is threefold:

• Qualify the SPH-FE coupling strategy, especially in terms of energy and then develop some schemes ensuring a good compromise among stability, accuracy and computation time.

• Quantify the influence of the different physical phenomena involved in aquaplaning problems.

• Adapt the SPH models to perform simulations of the complete problem.

Thesis organization

This manuscript is organized into several parts: first, the SPH method is introduced from the basics to advanced techniques, using the background of the LHEEA laboratory. The main advantages and limitations of the method are discussed. Then, a general introduction of the FE method is provided, before introducing the coupling strategy between SPH and FE methods. The way to impose the fluid loading and the deformable body treatment is detailed for different coupling algorithms. After these preliminary chapters stating the presentation of the numerical methods and the coupling strategy, an energetic study is conducted. The SPH equations are written focusing on energy conservation for the discrete particle system, especially on the interface coupling energy. This part aims at clarifying how the energy is dissipated/created at the fluid-structure interface, and at comparing this "artificial" energy with other energy components. Then, the investigations focus on computational cost improvements. The fluid-structure coupling strategy is validated on different 3D test cases. A limited number of studies is present in the literature for this specific topic. Three different experiments have been reproduced in this work.

The second topic to be thoroughly investigated is the general problem of a rolling tire on a puddle of water. A model is proposed to represent the aquaplaning phenomenon on a perfectly smooth ground.

Chapter 2

Smoothed Particle Hydrodynamics method

This chapter is dedicated to the presentation of the Smoothed Particle Hydrodynamics (SPH) method, especially through the state of the art in the SPH community, and more particularly at the LHEEA Laboratory. The SPH method was initially proposed by [START_REF] Lucy | A numerical approach to the testing of the fission hypothesis[END_REF] and [START_REF] Gingold | Smoothed Particle Hydrodynamics: theory and application to non-spherical stars[END_REF] at the end of the 70's in astrophysics.

Then [START_REF] Monaghan | Simulating free-surface flows with sph[END_REF] proposed an adaptation to free surface flow simulations. The range of applications of this method is large: fluid and solid mechanics, violent impacts, free surface flows, complex interfaces, multi-physics, etc.

The main characteristic of the method resides in its meshless approach and on an explicit solution. The scheme is therefore relatively easy to implement. It is based on an estimation of the discrete values of pressure and velocity derivatives from a set of particles at one specific instant. Then an explicit time integration is used to update positions, velocities, densities and pressures to the next instant, while avoid the solution of a Poisson equation. This feature provides a great flexibility and robustness to the method, from physical as well as algorithmic point of views.

Navier-Stokes Equations

Governing equations

Several approaches can be adopted to model incompressible flows in the field of computational fluid dynamics. It is possible to neglect the compressible effects and to restrain the compressible Navier-Stokes equations to model incompressible flows, i.e. div v = 0. Another approach to model an incompressible flow situation is to consider the fluid as weakly-compressible. In a Lagrangian description the Navier-Stokes equations can be expressed in the domain Ω by:

∀ x ∈ Ω          d x dt = v, dρ dt = -ρdiv( v), d v dt = g + divσ ρ , p = f (ρ), (2.1)
CHAPTER 2. SMOOTHED PARTICLE HYDRODYNAMICS METHOD where x and v are respectively the position and velocity vectors, ρ is the density and p is the pressure. g represents the external forces vector, here the gravity only. The last line of the set of equations is a state equation f , which thermodynamically closes the system of equations. In this work, it links the pressure and the density exclusively. Finally, σ is the stress tensor for a Newtonian fluid:

σ ij = (-p + λd kk )δ ij + 2µd ij , (2.2)
d ij = 1 2 ( ∂ v i ∂x j + ∂ v j ∂x i ), (2.3)
λ is the bulk viscosity, µ the kinematic viscosity and d the strain rate tensor. The conservative equations can be expressed in the following compact form:

∂ φ ∂t + div( F c -F υ ) = S. (2.4)
where φ is the vector of conservative variables, F c and F υ stand respectively for conservative and viscous flux tensors respectively and S is a source term:

φ = ρ ρ v , F c = ρ v ρ v ⊗ v + pI , (2.5)
F υ = 0 τ , S = 0 ρ g . (2.6)
where τ is the viscous stress tensor:

τ ij = λd kk δ ij + 2µd ij . (2.7)
Following the hypothesis of a constant dynamic viscosity ν = µ ρ and a quasi-incompressible fluid, the viscous term can be simplified and the Navier-Stokes equations become:

∀ x ∈ Ω          d x dt = v, dρ dt = -ρdiv( v), d v dt = g -∇(p) ρ + ν ∆( v), p = f (ρ), (2.8)

Equation of state

In order to solve the previous system Eq. (2.8) an equation of state is used to link pressure and density. In the literature a large variety of laws exist depending on the targeted applications. In this work the Tait equation of state is used (TAIT, 1888):

p -p 0 = ρ 0 C 2 0 γ κ ρ ρ 0 γκ -1 , (2.9)
where ρ 0 and C 0 are respectively the nominal density and sound speed of the fluid and p 0 is the reference pressure. C 0 controls the stiffness of the equation, and therefore the compressibility of the flow. The atmospheric pressure is often considered as, p 0 = 0. γ κ is the polytropric constant of the fluid, for water γ κ = 7 and for air γ κ = 1.4. This equation is based on the Laplace equation p ρ γκ = cst for adiabatic transformations of an ideal gas. The use of the Tait equation of state avoids solving the energetic balance in the system.

SPH FORMALISM

ALE formulation of Navier-Stokes equations

Previous equations were presented using a Lagrangian point of view where the temporal evolution of the fluid is based on its initial position. Nevertheless they can also be expressed using an Arbitrary Lagrangian Eulerian (ALE) description, which generalizes the previous equations to any transport field v 0 . The ALE conservative equations can be expressed as:

L v 0 (φ) + div( F c -F υ -v 0 ⊗ φ) = S, (2.10)
where L v 0 is the transport operator defined as: 2.11) which leads to: .12) where ω is the fluid volume.

L v 0 (φ) = ∂ φ ∂t + div( v 0 ⊗ φ), (
dω φ dt + ω div( F c -F υ -v 0 ⊗ φ) = ω S. ( 2
As a consequence, the system of equations can be rewritten as:

∀ x ∈ Ω        d x dt = v 0 , dρ dt = -ρdiv( v) -( v -v 0 ). ∇ρ, d v dt = g -1 ρ ∇. [ρ v ⊗ ( v -v 0 )] -1 ρ ∇(p) + ν ∆( v) + v ρ div(ρ v -ρ v 0 ), p = f (ρ).
(2.13)

The transport velocity can be equal to the flow velocity, leading to a Lagrangian description. Conversely, a transport velocity equal to zero corresponds to a Eulerian description.

SPH formalism

In this section the basic formalism of the SPH method is presented. The method has the main feature of being meshless. There are no pre-established connectivities among the interpolation points. The neighbourhood of each point evolves through time. The method is then particularly suitable for problems in which large deformations of the domain occur, possibly in the presence of complex interfaces and/or strong convection problems. Huge variations of the domain do not lead to particular problems unlike for mesh-based methods where mesh management procedures would be needed, leading to potentially strong computational costs. In return the SPH method usually suffers from low accuracy and is computationally demanding. This last issue can be however partially answered by massive parallelization thanks to the explicit nature of the method.

In a first approach only inviscid flows are considered, viscous flow considerations will be presented in a following section. This section is dedicated to introduce the SPH formalism for solving the Euler equations. The fluid is described by a set of particles (elementary fluid volumes), carrying the discrete field values within the fluid domain. This means that each particle represents a volume of fluid with its own mass, density, pressure and velocity. Time derivatives are then explicitly integrated in time.

CHAPTER 2. SMOOTHED PARTICLE HYDRODYNAMICS METHOD

Particle approach

In the SPH formalism, the space differential operators in the system Eq. ( 2.13) are approximated using a convolution product based on a weighting function. This is the fundamental difference with mesh-based methods.

Interpolation of a function φ

The approximation of a function φ at a location x can be expressed as: A radial kernel is used, so that:

< φ( x) >= D φ( y)W ( x -y)dV y , ( 2
W ( x -y, R) = W ( x -y , R). (2.15)
The compact property of the kernel support makes it equal to zero on the boundary ∂D as:

∀ y ∈ ∂D W ( x -y, R) = 0. (2.16)
The kernel function W converges towards the Dirac function as the size of its support tends to 0:

lim R->0 W ( x -y, R) = δ( x -y),
(2.17)

SPH FORMALISM

and W should satisfy:

D W ( x -y, R)dV y = 1. (2.18)
The gradient of the kernel ∇W presents an anti-symmetric property as:

∇ y W ( x -y, R) = -∇ y W ( y -x, R). (2.19)
Note finally that the approximation Eq. ( 2.14) is second order accurate as:

< φ( x) >= φ( x) + O(R 2 ).
(2.20)

Gradient interpolation of a function φ

The gradient convolution is defined as:

< ∇ x φ( x) >= D ∇ y φ( y)W ( x -y, R)dV y . (2.21)
Using an integration by parts and the anti-symmetric property of ∇W , it becomes:

< ∇ x φ( x) > = D ∇ y φ( y)W ( x -y, R)dV y = ∂D φ( y)W ( x -y, R) ndS - D φ( y) ∇ y W ( x -y, R)dV y = - D φ( y) ∇ y W ( x -y, R)dV y = + D φ( y) ∇ x W ( x -y, R)dV y . (2.22)
Approximation Eq. ( 2.22) is valid exclusively when the kernel support is complete since in this case:

∂D φ( y)W ( x -y, R) ndS = 0. (2.23)
This property is lost when the kernel support is truncated (close to the boundary of the domain for instance), leading to the need for taking this surfacic term into account.

Kernel functions

The choice of the kernel function is crucial. We introduce here two different kernel functions, but note that a wide variety can be found in the literature.

CHAPTER 2. SMOOTHED PARTICLE HYDRODYNAMICS METHOD

Cubic spline kernel

This kernel has been proposed by [START_REF] Monaghan | Simulating free-surface flows with sph[END_REF], and is often used by SPH practitioners.

W (q, R) = β    2 6 -2q 2 + 2q 3 if 0 ≤ q < 1 2 4 6 (1 -q) 3 if 1 2 ≤ q < 1 0 if q > 1 , (2.24)
and

∇W ( x, R) = β x x R    -4q + 6q 2 if 0 ≤ q < 1 2 -4(1 -q) 2 if 1 2 ≤ q < 1 0 if q > 1 , (2.25)
where q = x R and β is a normalization constant depending on the dimension of the problem, for 2D problems β = 120 7πR 2 and for 3D problems β = 24 πR 3 .

Wendland kernel

This kernel has been proposed by [START_REF] Wendland | Piecewise polynomial, positive definite and compactly supported radial function of minimal degree[END_REF].

W (q, R) = β 1 2 (1 -q) 4 (1 + 4q) if 0 ≤ q ≤ 1 0 if q > 1 , (2.26)
and

∇W ( x, R) = β x x R -10q(1 -q) 3 if 0 ≤ q ≤ 1 0 if q > 1 , (2.27)
where for 2D problems β = 14 πR 2 and for 3D problems β = 21 πR 3 . This kernel function shows very good stability and convergence properties compared to other kernel functions, [START_REF] Dehnen | Improving convergence in smoothed particle hydrodynamics simulations without pairing instability[END_REF]. It will be used for all the calculations in this thesis.

Discrete interpolation

The continuous formulation of spatial operators has been detailed. A discrete approach is then needed. Using the SPH formalism, the continuous domain is discretized using particles. An elementary volume ω i is provided to each particle i. Once the space of discretization introduced, it is possible to write the following discrete quadrature formula: (2.28) where N represents the number of particles in the support kernel D. Similar expressions can be defined for the discrete divergence and gradient operators:

D φ( y)W ( x -y, R)dV y ≈ N j=1 φ( x j )W ( x i -x j , R)ω j ,
< ∇ x i φ( x i ) >= N j=1 φ( x j ) ∇ x i W ( x i -x j , R)ω j .
(2.29)

Convergence and accuracy

Nevertheless these interpolations (2.28) and ( 2.29) have poor accuracy orders. Indeed, the divergence and gradient operators of constant functions are not correctly estimated. The discrete operators must also ensure the momentum and mass conservation equations to be consistent with the physical solution. In this work, the formulation of [START_REF] Bonet | Variational and momentum preservation aspects of smooth particle hydrodynamics[END_REF] obtained by applying the Principle of Virtual Work (PVW) (also adopted by COLAGROSSI ( 2001)) is used:

< ∇ x v( x) > = < ∇ x v( x) > -v( x). < ∇ x 1 > ≈ N j=1 [ v( x j ) -v( x i )] . ∇ x i W ( x i -x j , R)ω j , (2.30) < ∇ x p( x) > = < ∇ x p( x) > +p( x) < ∇ x 1 > ≈ N j=1 [p( x j ) + p( x i )] ∇ x i W ( x i -x j , R)ω j . (2.31)
For the rest of this thesis ∇ x i W ( x i -x j , R) and p( x i ) will be respectively noted ∇ i W ij and p i .

However the discrete approximations Eq. ( 2.28) and Eq. ( 2.31) do not fulfill exactly the previous integral properties because:

N j=1 W ij ω j = 1, (2.32) N j=1 ∇ i W ij ω j = 0. (2.33)
Theoretically an infinite number of particles would be required in the kernel support so that:

R ∆x → ∞, (2.34)
where ∆x is the characteristic distance between two particles. In practice this condition is never verified, and a constant ratio R ∆x is usually used in the range 2 < R ∆x < 4. Theses wrong properties are limited in the case of a Cartesian point distribution. However the Lagrangian nature of the SPH method generally leads to arbitrary distribution of points (see Fig. 2

.2).

For the integral interpolation to converge, an additional convergence property is also required: R → 0, (2.35) as R is a direct measure of the interpolation. This double convergence property cannot be respected in practice in the general case, so that various corrective procedure have been designed. For more informations about the accuracy and the precision of the SPH interpolation, please refer for instance to [START_REF] Quinlan | Tuncation errors in mesh-free particle methods[END_REF].

Shepard correction

The Shepard correction (SHEPARD ( 1968)) is useful to maintain a 0 th accuracy order for the discrete approximation. This corrective expression is:

< φ i >= N j=1 φ j W S ij ω j , (2.36)
where

W S ij = W ij N k=1 W ik ω k .
(2.37)

Renormalization

This second correction introduced by [START_REF] Randle | Smoothed particle hydrodynamics: Some recent improvements and applications[END_REF] proposes to increase the order of the gradient operator approximation to ensure the following discrete property:

N j=1 ( x j -x i ) ⊗ ∇ RL i W ij ω j = I d , (2.38) with ∇ RL i W ij = B i ∇ i W ij .
(2.39)

A renormalization matrix B i for each particle i is introduced:

B i = N j=1 ( x j -x i ) ⊗ ∇ i W ij ω j -1
.

(2.40)

WEAKLY-COMPRESSIBLE SPH

The gradient is then approximated using: (2.41) This technique leads to exact interpolation of constant and linear fields. Since the symmetric property is lost in Eq. ( 2.41), a symmetric renormalization can also be obtained by considering: .42) which leads to the following gradient operator:

< ∇ i φ i >= B i N j=1 [φ j -φ i ] ∇ i W ij ω j .
B ij = B i + B j 2 , ( 2
< ∇ i φ i >= N j=1 [φ j -φ i ] B ij ∇ i W ij ω j .
(2.43)

Weakly-Compressible SPH

Now the system (2.8) can be discretized using the discrete divergence and gradient operators defined in the last section. The concept of conservation is important to build a SPH scheme as stable as possible. Here, a weakly compressible approach is chosen, known as Weakly-Compressible SPH (WCSPH).

The operators previously defined are mainly centered in space VILA (1999). The explicit solution adopted therefore leads to an unconditionally unstable scheme. A stabilization method is necessary to guaranty the convergence of this system. In order to obtain a stable method, a numerical dissipation is introduced into the scheme. The main subtlety resides in the type of numerical dissipation to introduce and how to regulate the amount of energy to dissipate. Especially, the solution adopted should be consistent. Different techniques exist to design this dissipation. In this section, the discrete equations are provided for three formulations targeting inviscid fluids (Euler equations).

Scheme base on artificial viscosity

The first stabilizing technique was introduced by MONAGHAN (1992). A numerical dissipation is added in the momentum conservation equation through the Von-Neumann Richtmyer artificial viscosity, Π ij , adapted to the SPH formalism.

ρ i d v i dt = ρ i g - N j=1 [p j + p i + ρ j ρ i Π ij ] ∇ i W ij ω j , (2.44)
where

Π ij = -α AV R ij C ij ρ ij x ij . v ij x ij 2 + R 2 ij if v ij . x ij < 0 Π ij = 0 otherwise , (2.45) CHAPTER 2. SMOOTHED PARTICLE HYDRODYNAMICS METHOD where C i is the sound speed of the particle i. R ij = R i +R j 2 , C ij = C i +C j 2 , ρ ij = ρ i +ρ j 2 , v ij = v j -v i and x ij = x j -x i .
The coefficient α AV is test case dependent as it depends on the dynamics and the nature of the flow studied. Classically the value is taken in the range 0.01 < α AV < 0.1. The parameter = 0.01 avoids the division by zero when two particles coincide.

This stabilization technique is simple, easy to implement, and does not add prohibitive computational costs. This term brings the desired upwinding to the solution (depending on α AV ). These reasons make the use of this technique popular in the SPH community. Nevertheless, this technique quickly reaches its limits. Indeed, the coefficient α AV must be chosen large enough to avoid instabilities in the solution but should be as close as possible to zero to limit the total amount of numerical diffusion. Furthermore the amount of diffusion needed for stabilization purpose is not necessarily the same everywhere in the flow. A large dissipation is therefore often introduced without being actually needed. The α AV should be tuned for each case and may provide noisy pressure fields. Hence other techniques have been developed in order to improve the results, especially for pressure field regularity purposes.

δ-SPH scheme

This scheme was first proposed by MOLTENI andCOLAGROSSI (2009), and[START_REF] Ferrari | A new 3D parallel SPH scheme for free-surface flows[END_REF]. The present formalism was introduced by ANTUONO et al. ( 2010), and then analysed in [START_REF] Antuono | Numerical diffusive terms in weakly-compressible SPH scheme[END_REF][START_REF] Antuono | Energy balance in the δ-SPH scheme[END_REF]. The main idea is to add a diffusive term in the mass conservation equation. This term helps in increasing the pressure field regularity by reducing spurious high-frequency noises.

∀ x ∈ Ω              d x i dt = v i , dρ i dt = -ρ i N j=1 [ v j -v i ] ∇ i W ij ω j + N j=1 ψ ij ∇ i W ij ω j , ρ i d v i dt = ρ i g -N j=1 [p j + p i -Π ij ] ∇ i W ij ω j , p i f (ρ i ), (2.46)
where

Π ij = α δ RC 0 ρ 0 2 x ij . v ij x ij 2 , (2.47)
and

ψ ij = δ δ RC 0 (ρ j -ρ i ) - 1 2 ( ∇ RL ρ i + ∇ RL ρ j ). x ij x ij x ij 2 , (2.48)
δ δ and α δ control the amount of diffusion to introduce in the mass and momentum equations respectively, usually δ δ = 0.1 and α δ = 0.01. This variant is similar to the Monaghan formulation. Nevertheless it is more robust and provides much more regular pressure fields in a large range of cases. It also preserves the mass, the linear and angular momenta, with reasonable computational costs.

Riemann-SPH scheme

This last formulation is the one actually used throughout this study. It is based on the solution of Riemann problems at each interface separating a pair of particles. This technique was firstly introduced by Vila VILA (1999), and then adopted by [START_REF] Guilcher | Contribution au développement d'une méthode SPH pour la simulation numérique des interactions houle structure[END_REF], [START_REF] Ferrari | A new 3D parallel SPH scheme for free-surface flows[END_REF][START_REF] Quinlan | Development of the meshless finite volume particle method with exact and efficient calculation of interparticle are[END_REF]) and OGER et al. (2016b). It is based on the conservative Euler equations Eq. ( 2.4) with a formalism close to the Finite Volume Method for compressible flows. Following VILA (1999), the conservative form of the governing equations Eq. ( 2.4) for an inviscid flow (i.e. F υ (φ) = 0) can be discretized as:

∀ x ∈ Ω        dx i dt = v i , dω i dt = ω i (div v) i , d(ω i φ i ) dt + ω i ∇ i (F c (φ)) = ω i S i , p i = f (ρ i ).
(2.49)

Then using ∇F c (φ) = ∇F c (φ) + F c (φ) ∇1, (2.50) d(ω i φ i ) dt + N j=1 ω i ( F c i + F c j ). ∇ i W ij ω j = ω i S i . (2.51)
In this approach no artificial diffusion is directly added in the mass or momentum equations. The main idea is to solve a Riemann problem between each pair of particles. The problem is solved at mid-distance between the particles i and j with an interface velocity

v ij = v i + v j 2
and a normal direction

n ij = ∇ i W ij ∇ i W ij
. The left and right states of the Riemann problem correspond respectively to φ i and φ j . The centered flux ( F c i + F c j ) is then replaced by twice the Riemann problem solution F E :

dω i φ i dt + N j=1 ω i 2 F E ( φ i , φ j , n ij , v i , v j ). ∇ i W ij ω j = ω i S i . (2.52)
Finally the set of governing equations writes: (2.53) where ρ e ,P e and v e are the Riemann problem solutions (exact or linearized Riemann problem). In this scheme the particle mass is not exactly conserved; a symmetric mass exchange is observed between particles i and j since:

∀ x ∈ Ω              dx i dt = v i , dω i dt = ω i N j=1 2( v ij -v i ). ∇ i W ij ω j , dm i dt = -N j=1 2ρ e ( v e -v ij ). ∇ i W ij ω i ω j , d(m i v i ) dt = ω i ρ i g -N j=1 2[ρ e v e ⊗ ( v e -v ij ) + P e I] ∇ i W ij ω j ω i , p i = f (ρ i ),
-2ρ e ( v e -v ij ). ∇ i W ij ω i ω j = +2ρ e ( v e -v ij ). ∇ j W ji ω i ω j , (2.54)
ensuring the global mass conservation of the system.

CHAPTER 2. SMOOTHED PARTICLE HYDRODYNAMICS METHOD

In the present work, we preferred a variant of this formulation proposed by [START_REF] Leduc | Improvement of multiphase model using preconditioned riemann solvers[END_REF], where the local mass conservation is ensured. It implies a purely Lagrangian description of the medium, i.e. dm i dt = 0. As a consequence, the interface velocity v ij is chosen equal to v e :

∀ x ∈ Ω              dx i dt = v i , dω i dt = ω i N j=1 2( v e -v i ). ∇ i W ij ω j , dm i dt = 0, dm i v i dt = ω i ρ i g -N j=1 2P e ∇ i W ij ω j ω i , p i = f (ρ i ), (2.55)
This type of approach is closer to the actual physics of the flows. The artificial viscosity is replaced with the numerical diffusion embedded within the Riemann solution. It naturally depends on the flows according to the acoustic propagation. Furthermore no coefficient needs to be tuned to regulate the amount of dissipation. However it demands larger computational resources than the previous schemes.

MUSCL

The Monotonic Upstream-centered Scheme for Conservative Laws (MUSCL) proposed by [START_REF] Van Leer | Towards the ultimate conservative difference scheme. V a second order to godunov method[END_REF] is used to reconstruct the variables φ linearly at the i-j interface for accuracy issues. The physical solution is approximated using linear piecewise functions instead of constant piecewise functions. It is based on local estimation of the gradient function using the following renormalized operator B defined in Section 2.2.4.2.

< ∇ i φ i >= N j=1 (φ j -φ i )B i . ∇ i W ij ω j .
(2.56)

Note that the non-symetric renormalization is preferred here for accuracy purposes. The left and right states of the Riemann problem, φ L ij and φ R ji , are defined then as: .58) where

φ L ij = φ i + ( x ij -x i ). ∇ i φ i , (2.57) φ R ji = φ j + ( x ij -x j ). ∇ j φ j . ( 2
x ij = x i + x j
2 . However the above gradients cannot be directly used for monotonicity preservation purposes. These gradients are therefore limited (as proposed in the original MUSCL scheme) using slope limiters to obtain the left and right states of the Riemann problem φ l ij . Here the minmod limiter is used as:

φ l ij = φ i + β l ij (φ j -φ i ), 0 < β l ij < 1, (2.59)
where

β l ij = minmod(1, ∇ j φ j .( x ij -x i ) φ j -φ i ), minmod(a, b) = min(0, max(a, b)) + max(0, min(a, b)).
(2.60)

The use of this limiter strongly decreases the numerical oscillations in the solution.

TIME STEP AND TEMPORAL INTEGRATION

Time step and temporal integration

The weakly-compressible SPH formalism uses an explicit time integration. The main advantage resides in the fact that no linear system needs to be solved. Nevertheless such a scheme is known to be conditionally stable, so that the time step ∆t must respect some conditions.

Courant-Friedrich-Levy condition

The condition on the time step size is the Courant-Friedrich-Levy (CFL) condition. This conditions needs to be verified for all particles, and can be expressed as:

∆t i ≤ k CF L min i∈D R i C i , (2.61)
where C i is the sound speed of particle i and k CF L is the Courant number, chosen here as 0.375. In an isentropic transformation the sound speed C is defined as:

C = ∂p ∂ρ . (2.62)
For the Tait equation of state it becomes:

C = γ κ p ρ . (2.63)
If viscosity is considered, the time step needs to respect additionally another condition based on the viscous diffusion:

∆t ≤ k υ CF L min i∈D R 2 ρ i µ , (2.64)
where k υ CF L is the diffusion number, heuristically chosen as 0.125. This condition is dominant for low Reynolds numbers (Re ≈ 1).

Sound speed and weakly-compressible flows

A particular feature of the SPH method is to use a compressible model to approach an incompressible fluid flow. Compressible effects are supposed to be negligible when the Mach number Ma is smaller than 0.1. The nominal sound speed C 0 is therefore chosen to ensure small density variations. Note that it does not correspond to a strictly incompressible flow (an incompressible model corresponds to exactly null density variations). Nevertheless, it ensures that the compressible effects remain in the acoustic domain.

Eq. ( 2.61) emphasizes the need for choosing C 0 as low as possible while keeping Ma< 0.1, in order to maximize the time step size and to limit the CPU costs. The maximal time step is inversely proportional to the maximal sound speed. It is then needed to find a compromise between computational cost and incompressible solution.

CHAPTER 2. SMOOTHED PARTICLE HYDRODYNAMICS METHOD 2.4.3 4 th order Runge-Kutta scheme A 4 th order Runge-Kutta scheme is used in the present work. The scheme is composed of four steps:

φ n+1 = φ n + 1 6 ∆t(k 1 + 2k 2 + 2k 3 + k 4 ) (2.65)
with the four increments k 1 , k 2 , k 3 , k 4 expressed as:

k 1 = f RK (t n , φ n ) (2.66) k 2 = f RK (t n + 1 2 ∆t, φ n + 1 2 k 1 ∆t) (2.67) k 3 = f RK (t n + 1 2 ∆t, φ n + 1 2 k 2 ∆t) (2.68) k 4 = f RK (t n + ∆t, φ n + ∆tk 3 ) (2.69)
where φ n+1 et φ n represents the variables at time t n+1 and t n respectively and ∂φ ∂t = f RK (t, φ).

Boundary conditions

The problems encountered in fluid mechanics usually require the imposition of boundary conditions. As presented in Fig. 2.3, several kinds of conditions can be encountered:

• free surface

• solid walls • inflow/outflow boundaries ∂Ω = ∂Ω F S ∪ ∂Ω W ALL ∪ ∂Ω OU T ∪ ∂Ω IN .
(2.70) The treatment of boundary conditions is clearly one of the most challenging part of the SPH method and is still a very active research topic, [START_REF] Ferrand | Unified semi-analytical wall boundary conditions for inviscid, laminar or turbulent flows in the meshless sph method[END_REF][START_REF] Mayrhofer | Unified semi-analytical wall boundary conditions in sph: analytical extension to 3-d[END_REF]. Near a boundary, the kernel support is truncated and needs to be treated carefully.

BOUNDARY CONDITIONS

Wall conditions

Two different wall conditions can be applied, depending on the physical problem to solve and if viscous effects can be neglected or not. In the first case, a free-slip condition is considered on the boundary:

v f luid . n = v wall . n, (2.71)
with n the normal vector to the wall. If viscous effects are present, a no-slip condition must be considered by imposing the wall velocity equal to the adjacent fluid.

v f luid = v wall ,
(2.72) on the solid boundary.

Originally, the first technique proposed was based on Lennard-Jones forces between fluid particles and the wall (MONAGHAN, 1994) to complete the truncated part of the kernel support. In this manuscript, wall boundary conditions are imposed by using a ghost particle method and/or a normal flux method. Each of these methods has its own specificities, as presented hereafter.

Ghost particle method

The ghost-particle method has been introduced by LIBERSKY et al. (1993) and COLA-GROSSI and LANDRINI (2003a). This technique is similar to the ghost-cell method used in mesh-based methods. The main idea is to complete the kernel support by creating ghost particles on the opposite side of the wall. At each time step, a mirror particle is created for each near-boundary particle: each particle having a truncated kernel creates a ghost particle, Fig. 2.4.

It is possible to express the discrete approximation of the function φ and its gradient as:

< φ i >= N j=1 φ j W ij ω j + Ng j=1 φ j W ij ω j ,
(2.73)

< ∇ i φ i >= N j=1 (φ j -φ i ) ∇ i W ij ω j + Ng j=1 (φ j -φ i ) ∇ i W ij ω j , (2.74)
where N g is the number of ghost particles in the kernel support. These new particles take their physical variables from their mother particles: (2.76) in absence of gravity and where G(i) is the ghost particle of particle i. We refer the reader to DORING (2005) for more details in presence of gravity.

P G(i) = P i , (2.75) ρ G(i) = ρ i ,
A velocity is imposed to ensure the non-penetration of fluid particles through the wall as in Eq. ( 2.77) for free-slip conditions:

v G(i) = v i + 2( v wall . n -v i . n) n.
(2.77) For no-slip condition it becomes:

v G(i) = -v i + 2 v wall .
(2.78)

One of the huge advantages of this method is to not strongly affect the algorithmic structure. Once the ghost particles have been created, interactions are performed in the same way as for other particles. The main drawback of this method resides in the creation of the ghost particle itself. Indeed, for complex geometries (i.e. geometry with sharp angles for instance) the creation of ghost particles results in possibly large algorithmic complexity. The generalization of ghost particle creation procedures is indeed not straightforward at all. Then the second problem is to define properly the conservative variables of each ghost particle to ensure proper boundary conditions for complex boundaries (i.e in a corner). Furthermore all these operations can be computationally expensive. For more informations about ghost particle creation with singular geometries, please refer to [START_REF] Oger | Aspects théorique de la méthode SPH et applications à l'hydrodynamique à surface libre[END_REF]. A variant of this method involving fixed ghost particles (only calculated at the beginning of the simulation) was proposed by MARRONE et al. (2011).

Normal flux method

To overcome the limitations of ghost particles on complex geometries, an alternative method has been developed: the Normal Flux Method (NFM) MARONGIU ( 2007), [START_REF] Marongiu | Free surface flows simulations in Pelton turbines using an hybrid SPH-ALE method[END_REF][START_REF] De Leffe | Modélisation d'écoulements visqueux par méthode SPH en vue d'application à l'hydrodynamique navale[END_REF]. In this technique the kernel support of the particle is not completed, Fig. 2.5. Instead the Normal Flux Method introduces the calculation of the surfacic term in the SPH gradient approximation:

< ∇ x φ( x) >= D φ( y) ∇ x W ( x -y)dV y + ∂D φ( y)W ( x -y) ndS y .
(2.79) The formalism presented in DE LEFFE ( 2009) is adopted in the present work. It first applies a correction to the truncated kernel. The new approximation for a function φ becomes:

< φ( x) >= 1 γ( x) D φ( y)W ( x -y)dV y , (2.80)
where

γ( x) = D W ( x -y)dV y .
(2.81)

The Shepard correction term γ is introduced to ensure the first order accuracy. This correction is analytically estimated according to [START_REF] Leroy | Unified semianalytical wall boundary conditions applied to 2-d incompressible SPH[END_REF]), VIOLEAU et al. (2014). Using the kernel property, far from the boundary we should observe:

γ = 1.
(2.82)

Using the integration by parts, the gradient function therefore becomes:

< ∇ x φ( x) > = 1 γ( x) D ∇φ( y)W ( x -y)dV y = 1 γ( x) D φ( y) ∇ x W ( x -y)dV y + 1 γ( x) ∂D φ( y)W ( x -y) ndS y . (2.83) CHAPTER 2.

SMOOTHED PARTICLE HYDRODYNAMICS METHOD

Note that it can alternatively be written as:

< ∇ x φ( x) > = D φ( y) ∇ x W ( x -y)dV y + 1 -γ( x) γ( x) D φ( y) ∇ x W ( x -y)dV y + 1 γ( x) ∂D φ( y)W ( x -y) ndS y . (2.84)
where the first term corresponds to the "standard" formulation. The solid wall therefore introduces both surfacic and volumic contributions in the gradient calculation.

The discrete set of equations (2.55) using a Riemann-SPH approach becomes:

                                               d x i dt = v i , dρ i ω i dt = 0, dω i dt = 1 γ i N j=1 2ω i ( v e -v i ). ∇ i W ij ω j + 1 γ i Np j=1 2ω i ( v e -v i ). n j W ij S j , dω i ρ i v i dt = ω i ρ i g i -1 γ i N j=1 2P e ∇ i W ij ω j ω i -1 γ i Np j=1 2P e n j W ij S j ω i , p i = f (ρ i ), (2.85)
where n j and s j represent the normal and the area of the surface element j respectively. N p is the number of surface elements inside the kernel support. The flux at the wall is solved using a partial Riemann solver, as explained in DUBOIS (1982). The theoretical state at the boundary is estimated from the left state and the boundary conditions, which leads to find the density ρ e verifying v. n = 0. Then a MUSCL reconstruction is applied, but without a classical limiter according to DE LEFFE (2009). Nevertheless, some limitations based on the weakly-compressible assumption are imposed on the density and velocity variables in order to maintain the physical coherence:

   ρ l j = ρ j if 0.99ρ i < ρ j < 1.01ρ i ρ l j = ρ i otherwise , (2.86)    v l j = v j if v j < 0.1C i v l j = v i otherwise . (2.87)
The main advantage of this method is to naturally take into account complex geometries. It also gives very good results DE LEFFE (2009). The scheme is unchanged inside the fluid 2.5. BOUNDARY CONDITIONS domain (where γ = 1). Note that this scheme is no longer conservative, due to γ. Indeed there is no reciprocity for the interactions between two particles i and j. Note also that 3D simulations performed in the following chapters use γ = 1 in all the fluid domain, due to some difficulties to correctly estimate γ everywhere.

Free surface conditions

The presence of a free surface introduces the need for both kinematic and dynamic conditions. The kernel support is partially filled in the free surface area, and thus the interpolation loses its accuracy, as explained in Section 2.3. Moreover the position of the free surface and its local normal are unknown (and not easy to determine).

Free surface dynamic condition

The continuity of the stress vector must be ensured through the free surface. In monophasic inviscid simulations, this condition means imposing the atmospheric pressure p 0 at the surface. It seems difficult to impose such a condition on an unknown interface (see Fig. Nevertheless the dynamic condition is constantly applied (in an integral way), even if no terms appear directly in the calculation COLAGROSSI et al. ( 2009) (i.e. p 0 = 0). The advantage of the method lies in the fact that it is not necessary to know the position of the free surface for applying the dynamic condition. It is important to notice that this natural dynamic condition imposition is valid exclusively if no correction is applied to the kernel gradient. Indeed, the boundary condition is not satisfied when applying a Shepard or a renormalization correction. For those cases it is necessary to conserve a non corrected scheme in near free surface area. The implicit free surface dynamic condition gives very good results concerning the shape of the free surface and there is no limitation concerning free surface reconnection problems.

Free surface kinematic condition

This second condition for the free surface imposes the local velocity to be equal to the local fluid velocity. So that a particle at the free surface must stay on it. It is recognized that the Lagrangian approach of the SPH method naturally fulfils this condition, [START_REF] Colagrossi | Theoretical considerations on the free-surface role in the smoothed-particle-hydrodynamics model[END_REF]. Note that the free surface theoretically passes at the volume boundary of these particles, which means around a distance ∆x 2 from the particle location.

Considerations on viscous flows

Modelling viscous flows implies the need for adding the viscous flux tensor into the Euler equations. The most popular formulations used in the SPH method are detailed here. For more details about viscous flows using the SPH formalism, we refer the reader to DE LEFFE (2009).

Morris et al. formulation

This formulation was introduced by MORRIS et al. (1997). It is based on a combination of the SPH divergence operator and a finite difference scheme for the gradient:

∆ i v i = 2 N j=1 ( x j -x i ). ∇ i W ij x ij 2 ( v j -v i )ω j . (2.88)
This formula conserves the linear momentum. Nevertheless this operator presents some inconsistency issues in some configurations with free surface flows [START_REF] Colagrossi | Theoretical analysis and numerical verification of the consistency of viscous smoothedparticle-hydrodynamics formulations in simulating free-surface flows[END_REF].

Monaghan and Gingold formulation

This formulation proposed by [START_REF] Monaghan | Shock Simulation by the Particle Method SPH[END_REF] has a similar form to the one used for the artificial viscosity used to stabilize the system.

∆ i v i = K N j=1 ( x j -x i ).( v j -v i ) x ij 2 ∇ i W ij ω j , (2.89)
where K depends on the dimension of the investigated problem; for 1D: K = 6, for 2D: K = 8 and for 3D: K = 10. A relationship can be determined between the coefficient of artificial dissipation α AV and the dynamic viscosity: .90) This operator conserves the linear and angular momentum. Furthermore it does not suffer from inconsistency issues at the free surface. However, it is not possible to use this formulation with the Normal Flux Method (NFM) for wall boundary conditions [START_REF] De Leffe | Modélisation d'écoulements visqueux par méthode SPH en vue d'application à l'hydrodynamique navale[END_REF]. This formulation is retained in this work due to the presence of a free surface when the viscous effects are considered.

µ = α AV hCρ K . ( 2

PARTICLE SHIFTING

Particle shifting

The accuracy of the SPH operators is directly linked to the particle distribution inside the kernel support. The Lagrangian nature of the SPH method makes the particles to follow the fluid trajectories, which may affect the quality of the simulation. Particles may be very irregularly spaced and may present some strongly anisotropic distributions, penalizing the spatial interpolations. Paradoxically any gain in accuracy makes the particle to follow more precisely the Lagrangian flow trajectory, decreasing the spatial interpolation accuracy. Various methods exist (often called "particle shifting") in the literature for disordering particles and aiming at preventing the occurrence of these anisotropic particle structures. Two of them are presented here. Nevertheless these corrections cannot be applied everywhere in order to preserve the free surface kinematic condition. A free surface detection should therefore be performed. This is done by applying the method proposed by MARRONE et al. (2010).

X-SPH

The first particle shifting procedure was proposed by [START_REF] Monaghan | On the problem of penetration in particle method[END_REF] by slightly modifying the particle displacements to obtain a better distribution and therefore a better accuracy. The main idea relies on moving particles with their Lagrangian velocities but corrected by a small velocity perturbation, as:

d x i dt = v i + δ v i , (2.91)
with

δ v i = N j=1 m j ρ ij ( v j -v i )W ij , (2.92)
where

ρ ij = ρ i +ρ j 2 .
In practice this correction avoids particle clustering. is a parameter controlling the velocity variations, typically = 0.5. Nevertheless, the X-SP H method tends to affect the conservation and the consistency properties of the global scheme, due to the fact that this small velocity offset only addresses the particle positions without any consideration of the effects on mass and momentum equations, which is not exact from a mathematical point of view.

ALE-SPH

OGER et al. (2016b) proposed a specific transport velocity within the ALE formalism presented in Eq. (2.13) for particle shifting purpose. ∀ x ∈ Ω :

             dx i dt = v 0i , dω i dt = ω i N j=1 ( v 0j -v 0i ). ∇ i W ij ω j , dm i dt = -N j=1 2ρ e ( v e -v 0ij ). ∇ i W ij ω i ω j , dm i v i dt = ω i ρ i g -N j=1 2[ρ e v e ⊗ ( v e -v 0ij ) + P e ] ∇ i W ij ω j ω i , p i = f (ρ i ), , (2.93)
with v 0i = v i + δv i . (2.94)
The main issue in particle shifting procedures is to design the perturbation velocity vector δv i . OGER et al. (2016b) proposed a method based on the Mach number M a i of the particle.

δv i =      -U char i R i k if U char i R i k < 0.25 v i -0.25 v i k k
otherwise.

(2.95)

where k is defined as

k = N j=1
∇W ij ω j , and

U char = M a C 0 .
The ALE formalism allows imposing a shifting transport velocity field while preserving conservation and consistency properties of the scheme. Unlike most of the existing particle disordering/shifting methods, this formalism avoids the formation of anisotropic structures while presenting a full consistency with the original Euler or Navier-Stokes equations. This formalism therefore provides a gain in terms of accuracy, convergence and numerical diffusion.

Local refinement for the SPH method 2.8.1 Variable-h method

In the previous sections the kernel radius was considered as the same for all particles in the fluid domain. This can lead to expensive computational costs due to large numbers of particles for wide domains and fine spatial resolutions. The use of varying spatial resolutions can help overcoming these issues, by concentrating the accuracy in desired identified areas and degrading the resolution elsewhere. For instance, high accuracy is needed in an impact area, or any area where very local effects occur. The notion of varying spatial resolution in the SPH method was first been introduced by [START_REF] Hernquist | Treesph: a unification of SPH with the hierarchical tree method[END_REF], [START_REF] Nelson | Variable smoothing lengths and energy conservation in smoothed particle hydrodynamics[END_REF] by slowly varying the support radius R. This technique offers a great interest and approaches the local refinement techniques often used in mesh-based methods. It pursuits the kernel radius to be small in this desired area and bigger elsewhere. Nevertheless, this leads to the introduction of dR dt and ∇R terms in the discrete equations. To use this technique in a simplified way, the formalism presented by OGER et al. ( 2006) is adopted here. Each particle possesses its own kernel radius, which remains constant during all the simulation. It prevents from the introduction of dR dt terms ( dR dt = 0). Furthermore the kernel radius varies slowly in space in order to limit the magnitude of ∇R terms so as to neglect it. It limits the amount of errors generated through the spatial variation of R. Typically, the kernel radius is initialized with a variation lower than 3%. It stands as a good compromise between accuracy and computational costs. The discrete approximation of a function and its gradient is performed using a symmetric form by replacing R with

R ij = R i +R j 2
in ∇W ( x i -x j , R). This formulation conserves the reciprocity of the density variations and the linear and angular momentum. It also conserves the global properties of the scheme.

LOCAL REFINEMENT FOR THE SPH METHOD

Nevertheless, the ratio of 3% on the radius R can lead to large transition zone between the finest and coarsest areas. Furthermore, depending on the flow physics coarse particles can be mixed into the fine spatial resolution domain leading to severe numerical instabilities. To overcome these drawbacks, new methods based on local particle refinement have been developed.

Adaptive Particle Refinement

The Adaptive Particle Refinement method (APR) used in this work was initially proposed by [START_REF] Feldman | Dynamic refinement and boundary contact forces in SPH with applications in fluid flow problem[END_REF] and then improved by [START_REF] Barcarolo | Adaptive particle refinement and derefinement applied to Smoothed Particle Hydrodynamics method[END_REF] and [START_REF] Chiron | Couplages et améliorations de la méthode SPH pour traiter des écoulements multi-échelles temporelles et spatiales[END_REF]. The method proposed by CHIRON et al. ( 2017) is inspired by the Adaptive Mesh Refinement (AMR) of the mesh-based methods. A particle can be refined (parent particle) by being divided into a finite number of smaller particles (called child particles) following a uniform Cartesian pattern (see Fig. 2.7). The number of child particles of each parent particles is equal to 4 and 8 for 2D and 3D cases respectively. Two sets of particles are therefore introduced: c containing the child particles and p containing the parent particles, verifying:

∆x c = 0.5∆x p , (2.96)
R c = 0.5R p . (2.97)
A mass ratio λ c is also defined, providing the mass of each child particle:

m c = λ c m p , (2.98)
2D: λ c = 0.25, 3D λ c = 0.125 ensuring the mass conservation. In this approach, the two sets of particles do not interact directly with each other. The interpolations between these two sets are performed through the use of guard particles, similar to guard cells used in AMR. The refinement procedure obtained can be seen as the union of two sets of particles having each its own refinement level (see Fig. 2.8). At the beginning of each sub-step of the fluid temporal integration scheme, the guard particles estimate their conservative variables φ × using a Shepard interpolation:

• • • • • • • • • R p ∆x p • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • ∆x c R p R c
< φ × i >= N • j=1 φ j W ij N • k=1 W ik ω k ω j . (2.99)
where N • corresponds the SPH particles inside the support kernel.

Then the SPH child particles (represented by •) estimate their time derivatives using the classical SPH formalism on the child particles (represented by • and ×), as presented in Fig. 2.8. The same procedure is performed for the SPH parent particles (represented by •) using all parent particles (represented by • and ×). On the contrary, the variable fields of guard child particles (represented by ×) and guard parent particles (represented by ×) are estimated using the Shepard interpolation on the SPH particles (represented by • and •). The spatial operators are made between pairs of particles having compatible kernel radius. The kernel radius of guard parent particles located inside the refined area is adjusted to the child one.

In practice a box is defined inside the fluid sub-domain in which the parent particles have to be refined. A buffer zone of size R p is defined where guard child particles are created. Conversely, a SPH parent particles leaving the buffer zone to the refined area becomes a guard parent particle. This particle is maintained into the flow and does not interact with the flow. It uses the variable fields of the child particles to update its trajectory in the refined domain, and become again a SPH parent particle outside the refined domain. On the contrary, any child particle located out of the box is finally erased [START_REF] Barcarolo | Adaptive particle refinement and derefinement applied to Smoothed Particle Hydrodynamics method[END_REF].

The use of a buffer zone is necessary to prevent pressure discontinuities which may occur when a parent particle enters (or leaves) a refined domain, instantaneously increasing (or decreasing) the number of neighbour particles per particle. Nevertheless it is not sufficient to maintain the pressure field regularity at the interface when the dynamic flow is violent and the advective term is dominant. [START_REF] Chiron | Analysis and improvements of adaptive particle refinement (APR) through CPU time, accuracy and robustness considerations[END_REF] proposed to handle this problem by spatially reorganizing the particles using the particle shifting presented in Section 2.7. In the same article, they showed how this APR technique is computationally efficient and accurate.

Chapter 3

Finite element method

This chapter is dedicated to the presentation of the Finite Element (FE) method. For the last decades, numerical methods were mainly mesh-based. For structure modelling, FE is usually adopted [START_REF] Bathe | Finite Element procedures[END_REF]. This chapter gives a brief description of the method, for more details about FE theory, we refer the reader to [START_REF] Belytschko | Nonlinear Finite Elements for Continua and Structures[END_REF], [START_REF] Fish | First course in finite elements[END_REF].

The FE method is a well-known numerical tool for structure analysis. As for fluid mechanics, the solid mechanics is governed by differential equations for which analytical solutions do not exist. The FE formalism can be obtained using different approaches. FE solves differential equations characterizing the physics of the structure in a discrete way. The main advantage of this method resides in its maturity together with its capability to model complex geometric structures, materials, fractures, etc. A Lagrangian FE method is used here to agree with the SPH main feature. The Lagrangian description of motions is only valid if the structure undergoes large deformations without mesh distortion, which is expected in this work. In this thesis, three different sources of non-linearity are considered: non-linear geometry, non-linear material and contact problems. The first section provides the governing equations derived from the principle of virtual work for a general continuous solid and spatial discretization. The second section presents the discrete system in matrix form. Finally, the time integrators are discussed at the end of this chapter.

Weak form for updated Lagrangian formulation

Introduction

Let us consider the domain Ω s in which the equations read:

d 2 u dt 2 = g + divσ ρ s , (3.1) σ = A , (3.2)
where u is the solid displacement, ρ s is the density of the solid, ( u) is the strain tensor and A the matrix of the stress-strain constitutive relation.

CHAPTER 3. FINITE ELEMENT METHOD

Initial conditions in the solid domain Ω s are:

u(t = 0) = u 0 , (3.3) d u dt (t = 0) = v 0 . (3.4)
Different boundary conditions can be imposed on the solid boundary domain ∂Ω s :

∂Ω t ∪ ∂Ω v = ∂Ω s , (3.5) ∂Ω t ∩ ∂Ω v = ∅, (3.6)
where ∂Ω t and ∂Ω v are respectively the solid boundaries for the strain and the velocity constraints.

On

∂Ω v : v = v b , (3.7)
and on ∂Ω t :

σ. n b = T b , (3.8)
where n b and v b are respectively the normal and velocity vectors of the boundary and T b is the traction vector.

PVW approach

The continuous problem is presented through the virtual work principle, by multiplying equation (3.9) with a virtual displacement δ u before integrating over the solid domain.

Ωs δ u. d 2 u dt 2 - divσ ρ s -g dV = 0, (3.9)
and using the Green-Ostrogradsky theorem and the symmetricity σ. Note that this relation needs to be satisfied for any arbitrary δ u values. Thus n weighting functions δ u can be chosen to obtain a system of n equations. ∀i ∈ [1, n]:

Ωs δ u i .ρ s d 2 u dt 2 dV + Ωs σ : δ (δ u i )dV - ∂Ωt δ u i . T b dS - Ωs δ u i .ρ s gdV = 0, (3.14)
The choice of weighting functions is completely free provided they are linearly independent, and provides more or less accurate results.

Matrix formulation of the PVW

Once the PVW obtained, an approximate solution u * of u is estimated using a linear combination of n functions known as shape functions. A Galerkin approach can be used, the weight and shape functions are identical. The method then consists in weakening the previous integrals, satisfying for only n weighting functions. The approximation of u * can be expressed as:

u * = W q, (3.15)
where W is the matrix of shape functions and q the parameter vector of the approximation. Provided that the choice of shape functions and weighting functions leads to a well-conditioned system, its solution will provide us with an approximation of the problem. The weighting function is:

δ u = W δ q, (3.16)
Using a matrix representation of the tensors:

σ : δ = T σ (3.17) with T = < xx , yy , zz , 2 xy , 2 xz , 2 yz >, (3.18) σ T = < σ xx , σ yy , σ zz , σ xy , σ xz , σ yz > . (3.19)
It leads to:

T = L u * , (3.20) σ = A .
(3.21) CHAPTER 3. FINITE ELEMENT METHOD Assigning these terms to the variational formulation

δ q. Ωs W T ρ s W q + (LW ) T A(LW ) q dV = δ q. ∂Ωt W T b dS + Ωs W ρ s gdV . (3.22)
Thus the PVW reads:

M q + K q = F ext , (3.23)
where M and K are respectively the mass and stiffness matrices and where F ext is the generalized force vector.

M = Ωs W T ρ s W dV, (3.24) K = Ωs (LW ) T A(LW )dV, (3.25) F ext = Ωs W ρ s gdV + ∂Ωt W T b dS. (3.26)
The matrix formulation is discretized in the next section.

Finite element discretization

Spatial discretization

The main idea is to divide the solid domain Ω s into a finite number of elements.

Ω s = ∪ N F E j=1 Ω j , (3.27)
where Ω j is the volume of the element j. Contrary to the SPH method, all elements are connected by nodes. The approximate solutions are given at this mesh points. Due to the Lagrangian nature, nodes and elements move with the material, which is critical for history-dependent materials. The domain boundaries are given by the element edges. Thus, complex boundary shape unavoidably introduces some errors. The approximate solution is accurate if the space of functions is large enough to correctly represent the solution. The principle is based on:

1. the construction of the function δ u for each finite element.

2. the use of nodal variables as approximate parameters which allows to impose boundary conditions to moving problems.

FINITE ELEMENT DISCRETIZATION

Nodal approximation

The finite element method is based on the systematic construction of an approximation u * of u by subdomain. This approximation is built on the approximate values of the field at the nodes of the element. In the sub-domain Ω j for element j:

u * = N u n , (3.28)
where N is the matrix of shape functions of element j. Generally, shape functions are polynomial-based functions. There is the same number of functions than the number of interpolation points multiplied by the number of variables per point.

x n is the nodal variable associated to interpolation nodes of element j. Thus:

δ u = N δ u n , (3.29) T = LN u n , (3.30) σ = A = ALN u n . (3.31)
The PVW reads:

δ u n . M j ün + K j u n -F ext,j -F link,j = 0, (3.32)
Element mass matrix:

M j = Ω j N T ρ s N dV, (3.33)
Element stiffness matrix:

K j = Ω j (LN ) T A(LN )dV, (3.34)
External loads on the boundary:

F ext,j = Ω j N T {ρ s g} dV + ∂Ω t,j N T {T b } dS. (3.35)
Unknown efforts F link represent external mechanical actions to the element. The summation of liaison efforts is thus null. Note that for any arbitrary virtual displacement δ u Eq. (3.32) is satisfied leading to:

M j ün + K j u n = F ext,j + F link,j , (3.36)
Then by summing over all elements, the global mass and stiffness matrices can be determinated. Depending on the problem, different solvers can be used to compute the unknown displacements. A time integrator updates the different quantities and their derivatives according to an implicit or an explicit scheme.

Time integrator scheme

The choice of the time integration scheme is a key point in the Finite Element method (and especially for its coupling with SPH regarding stability). Several models can be used: implicit or explicit schemes. For fast dynamics, explicit schemes are classically preferred to implicit ones. In this work, we propose to test different time integration schemes in the structure solver to obtain the most stable and robust coupling strategy. Here we focus on two classical implicit schemes, Newmark and Hilber-Hughes-Taylor (HHT), and an explicit scheme using a central-difference integration rule.

Newmark implicit scheme

At time t n+1 equilibrium can be expressed at each mesh node as:

M d 2 u dt 2 n+1 + F n+1 int = F n+1 ext . (3.37)
Assuming that all variables are known at t n , the new equilibrium solution is obtained by using the Newmark scheme:

u n+1 = u n + ∆t d u dt n + ( 1 2 -β s )∆t 2 d 2 u dt 2 n + β s ∆t 2 d 2 u dt 2 n+1 , (3.38)
d u dt n+1 = d u dt n + (1 -γ s )∆t d 2 u dt 2 n + γ s ∆t d 2 u dt 2 n+1 .
(3.39)

The scheme stability depends on the values of the parameters β s and γ s . The scheme is unconditionally stable with

γ s ≥ 1 2 , (3.40) and β s ≥ 1 4 (γ s + 1 2 ) 2 .
(3.41)

A second order convergence is expected with

γ s = 1 2 , (3.42)
without numerical diffusion in the structure. Nevertheless, for

γ s > 1 2 , (3.43)
a numerical diffusion is introduced; it is usually needed to stabilized the solution if the structure is exposed to high frequency pressure loads. This is particularly relevant for coupling with the SPH method which usually imposes high frequency acoustic pressure noises due to its weakly-compressible nature. However, this diffusion can also affect the low frequency domain, thus harming the solution accuracy.

TIME INTEGRATOR SCHEME

Hilber-Hughes-Taylor scheme

To reduce high frequency oscillations without affecting the low frequency domain, the Hilber-Hughes-Taylor scheme has been proposed by [START_REF] Hilber | Improved numerical dissipation for time integration algorithms in structural dynamics[END_REF]. This approach introduces an offset within the equilibrium equation:

M d 2 u dt 2 n+1 + (1 + α s ) F n+1 int -α s F n int = (1 + α s ) F n+1 ext -α s F n ext , (3.44) with α s ∈] -1 3 , 0], related to γ s = 1 2 -α s , (3.45)
and

β s = (1 -α s ) 2 4 . (3.46)
This offset increases the numerical diffusion on the high frequency domain but reduces the numerical diffusion for low frequencies while ensuring a second order convergence. It can be noticed that for α s = 0, the HHT scheme corresponds exactly to the Newmark one with γ s = 1 2 and β s = 1 4 .

Central difference scheme

Equations of motion can be integrated in time using the following explicit central-difference integration rule.

d u dt

n+ 1 2 = d u dt n-1 2 + ∆t n+1 + ∆t n 2 d 2 u dt 2 n , (3.47)
u n+1 = d u dt n + ∆t n+1 d u dt n+ 1 2 .
(3.48)

It belongs to the family of Newmark scheme with β s = 0 and γ s = 0.5. The explicit integration rule is simple. Nevertheless, as for the fluid, central-difference operator is conditionally stable. The time step in the structure must respect the solid CFL condition which depends on material and mesh characteristics:

∆t < min L e C d , (3.49)
where L e is a characteristic length associated to an element and C d is the dilatational wave speed of the material. However, this wave speed is generally much higher than water speed of sound and it cannot be decreased without affecting the scheme stability. This leads to very small time steps compared to implicit schemes, thus increasing the computational time. Another problem lies in the need of numerical dissipation as shown in the following sections.

CHAPTER 3. FINITE ELEMENT METHOD

Discussion

This section has provided a brief introduction to the FE method. The method has not been presented in detail because our coupling is not based on hard coding within the FE model. Indeed, our aim is to develop a coupling strategy adaptable to any structure solver and problem. The FE solver is therefore considered here as a "black box", receiving fluid pressure loads and returning structure node positions and velocities. Nevertheless, the informations provided in this chapter are sufficient to understand the basics of the method and its limitations. For more information about the FE theory, please refer the related literature, such as BATHE (1995) and [START_REF] Belytschko | Nonlinear Finite Elements for Continua and Structures[END_REF]. The next chapter is dedicated to SPH-FE coupling strategy.

Chapter 4

SPH-FE coupling strategy for FSI problems

A coupling strategy between Smoothed Particle Hydrodynamics and Finite Element methods is proposed in this chapter, following the discussion in Section 1.2. It is dedicated to violent fluid-structure interaction modelling. The use of a Lagrangian meshless method for the fluid reduces the complexity of fluid-structure interface handling, especially in the presence of complex free surface flows such as the aquaplaning phenomenon. The coupling algorithms and information transfers are introduced here. Both convergence and robustness properties of the SPH-FE coupling are also discussed.

Coupling strategy

The discussion revolves around the coupling strategy adopted by [START_REF] Fourey | Développement d'une méthode de couplage fluide structure SPH Eléments Finis en vue de son application á l'hydrodynamique navale[END_REF][START_REF] Fourey | An efficient FSI coupling strategy between smoothed particle hydrodynamics and finite element methods[END_REF]. The present thesis focuses on SPH-FE coupling to model violent fluidstructure interaction problems involving complex free surface flows with deformable structures. The solid part is thus modelled with a classical FE method. The fluid part uses the SPH method to easily consider complex free surface flows. The meshless and Lagrangian natures of SPH naturally avoids the problem of fluid/solid grid compatibility. Neither contact algorithms between fluid and solid nor free surface tracking algorithms are needed. This property results in a significantly reduced complexity of the fluid-structure interface treatment. Another advantage of the SPH method resides in its compressible formulation, allowing for considering actual compressible effects when needed (cases involving high dynamics, fluid squashing...). Note also that this coupling strategy leads to a purely Lagrangian description in both media. Classical partitioned approaches are used, aiming at developing a coupling strategy that is robust and compatible with any kind of SPH and FE solvers. Here the deformable body behaviour is modelled using two generalist software: an open-source one Code Aster developed by EDF and a commercial one Abaqus developed by Simulia. The fluid part is modelled using SPH-flow, jointly developed by Ecole Centrale Nantes and NEXTFLOW Software (OGER et al., 2016a,b).

Of course, the choice of the solid time integration scheme is important. Explicit schemes are usually preferred to model very fast dynamics phenomena such as impacts. On the contrary, implicit schemes are usually used for less violent dynamics, where physical times are longer and require larger time steps. Note also that fast dynamics for structures happens usually with time scales which are much smaller than what can be considered as "fast dynamics" in fluids.

Besides, FOUREY (2012) showed that the solid time integration scheme can also be used to stabilize the coupled SPH-FE calculations. The weakly-compressible nature of the SPH method is responsible for acoustic pressure waves, which can affect the coupling stability. [START_REF] Fourey | Développement d'une méthode de couplage fluide structure SPH Eléments Finis en vue de son application á l'hydrodynamique navale[END_REF] showed that the introduction of sufficient dissipation on the high frequency domain maintains a stable coupling calculation, for example via the implicit FE time integration scheme. Thus, implicit and explicit schemes should respectively be used for the FE and SPH solver. The structural part is solved by using common implicit algorithms, such as the Newmark and HHT schemes, which introduce a sufficient dissipation on the high frequencies domain.

Deformable body considerations

Interaction which a deformable body differs from with a rigid body and requires therefore specific arrangements. This section aims at presenting the exchanges between SPH and FE solvers and explaining the interface treatment from each part.

Treatment of the deformable body

From the SPH algorithmic point of view, a deformable body at instant t is considered as a rigid wall, but with specific local deformation velocities (and node positions) provided by the FE solver. It is possible to choose any boundary method, ghost particle method or NFM (Fig. 4.1). No other treatment is prescribed to treat the deformable body interface. Especially, no contact algorithm dedicated to avoid material interpenetration is therefore needed, contrary to GROENENBOOM and CARTWRIGHT (2010) for instance. As proposed by [START_REF] Oger | Twodimensional SPH simulations of wedge water entries[END_REF], the pressure loads P k applied to each the FE structure mesh is computed using the SPH fluid pressure solution averaged from the near face boundary area.

Fluid loads on the deformable structure

More precisely, pressure loads are defined as an average of all particles seen by the wet body panel (Fig. 4.2) within a distance R from the panel, where R is the kernel radius:

P k = 1 N N i∈Aa P i . (4.1)
where A a is the averaging area.

The average pressure area is chosen proportional to R for spatial convergence purposes, since R is a convergence measure in the SPH method. As emphasized by [START_REF] Oger | Twodimensional SPH simulations of wedge water entries[END_REF], pressure averaging tends to smooth the irregularities near the wet body panel, which is important in the present work as pressure load irregularities may be responsible for coupling instabilities. A sufficient number of particles is necessary for each panel, of which the size has to be adapted in order to capture the physical load variations. The influence of the number of particles per panel onto the coupling robustness will not be directly treated here, but could be the purpose of a future work.

Coupling algorithms

Two different schemes are considered (see Fig. 4.3): the Conventional Parallel Staggered (CPS) procedure for parallel algorithm and the Conventional Sequential Staggered (CSS) for sequential algorithm. Both schemes are described by [START_REF] Farhat | Two efficient staggered algorithms for the serial and parallel solution of three dimensional nonlinear transient aeroelastic problems[END_REF].

The time step is the same for fluid and solid solvers. As the maximum fluid time step allowed is expected to be smaller than the solid one, SPH-Flow imposes its time steps to the FE solver. This condition is always satisfied using an explicit resolution in the fluid and an implicit one in the solid. In practice:

• for the CPS algorithm the fluid solver sends its time step value and pressure loads to the FE solver and receives the body node positions and velocities. Then both solvers simultaneously evolve from t n to t n+1 .

• for the CSS algorithm, fluid and solid solvers progress alternatively: the fluid solver receives the body node positions and velocities before progressing to t n+1 . Then, it sends its time step value and pressure loads to the FE solver and waits until the solid solver has reached t n+1 . The alternating operation of the CSS algorithm introduces a shift of the solution and improves the coupling stability. It also correspond to the algorithm that best fits with the 4.4. COUPLING ROBUSTNESS implicit nature of the solid resolution. Indeed the CPS algorithm makes the assumptions that the fluid loading F f luid/body does not vary much between two instants (i.e. F n+1 f luid/body ≈ F n f luid/body ). It is acceptable since the weakly compressible approach imposes very small time steps. Nevertheless, CPU costs are increased for the CSS scheme compared to the CPS one. Indeed, the calculation time of fluid and solid solvers are no longer hidden by a synchronous processing and are therefore added together. It is then important to establish which one of these two algorithms corresponds to the best compromise in terms of stability/CPU time.

CPS coupling algorithm

CSS coupling algorithm

Note that for both schemes no sub-iterations are performed to ensure velocity and pressure agreements at the fluid/solid interface. There is no convergence or stability control during the simulation, which differs from common practice with conventional mesh-based methods. Actually, the WCSPH method imposes very small time steps, so that only small variations of the different quantities are observed between two time steps. Thus, it naturally introduces a "natural" convergence in time which helps in stabilizing the coupling.

Recently LI et al. ( 2013) proposed a specific coupling strategy enforcing energy conservation at the interface. This technique is strongly intricated with the FE model, and therefore implies hard-coding within the solid solver. Conversely, the weak-coupling strategy retained in this study is flexible, easy to implement and compatible with any kind of SPH schemes and boundary condition treatment. The solid solver is considered as a "black box" receiving fluid loads and returning updated solid positions and velocities, so that any FE software can actually be used. No additional CPU time is introduced as no interface treatment is required.

Coupling robustness

A robustness analysis of this coupling strategy between SPH for the fluid and FE for the structure has been presented by [START_REF] Fourey | An efficient FSI coupling strategy between smoothed particle hydrodynamics and finite element methods[END_REF]. More particularly, the loss and gain in stability is studied according to various coupling parameters, and different coupling algorithms are considered. Validation analysis has been performed on 2D simulations of academic [START_REF] Scolan | Hydro-elastic behaviour of a conical shell impacting on a quiescentfree surface of an incompressible liquid[END_REF] and experimental [START_REF] Antoci | Numérical simulation of fluid structure interaction by SPH[END_REF] test cases.

The importance of adding some diffusion in the FE time integration scheme for maintaining the stability of the coupling has been emphasized. Acoustic waves tend to trigger the destabilization of the coupling. Note that acoustic waves occur for any fluid simulation solved using a weakly compressible approach (as in a real fluid). Instabilities in the structure solution lead to interface high frequency displacements of the interfacial nodes which are responsible for acoustic waves within the fluid. Subdividing the time step using finer spatial resolutions helps mitigating this phenomenon but only partially. Adding a small numerical diffusion in the FE time integration scheme (Newmark or HHT) tends to cancel this noise and allows stabilizing the coupling. We refer the reader to KRENK (2016) for more details about the dissipation inside the Newmark scheme. Note that a relatively strong damping (increasing with |α s |) does not introduce any phase shifting. Once the coupling strategy stabilized, it gives identical solutions whatever the value adopted for α s .

To highlight this sensitivity to high frequency fluid loading, a test case originally proposed by FOUREY ( 2012) is considered. It consists in a hydrostatic water column on an elastic plate according to the configuration presented in Fig. 4.4 and Tab. 4.1. The struc-CHAPTER 4. SPH-FE COUPLING STRATEGY FOR FSI PROBLEMS ture mesh is composed of 4 width-wise and 40 length-wise elements. The fluid is initially at rest on an aluminum plate which is built-in at its two extremities. This beam is initially not deformed since it is considered as not loaded yet with the fluid pressure forces (i.e. with the fluid weight).

We focus here more particularly on the midpoint displacement time history of the aluminum beam. The fluid-structure coupling has to find the equilibrium solution of this hydrostatic problem. Note that this case has also been treated by [START_REF] Li | A non-intrusive partitioned approach to couple smoothed particle hydrodynamics and finite element methods for transient fluid-structure interaction problems with large interface motion[END_REF] in order to show the stability property of their fluid-structure coupling algorithm compared to ordinary staggered methods, such as those used here (CPS and CSS). This test case is particularly interesting since the analytical solution of the problem is known, while analytical FSI validations are very difficult to find in the related literature. The structure is first initialized by performing a static calculation without any load except its own weight. At t = 0, the aluminum plate is suddenly uniformly loaded with the water column weight. The beam therefore starts to oscillate in response to this sudden change. These high frequency oscillations in the structure (period of about 4 ms) contrasts 4.4. COUPLING ROBUSTNESS with the acoustic wave frequency in the fluid, corresponding to 2H (wave propagation back and forth along the water height) at speed C 0 (period of about 80 ms). Two indicators are studied. The first is the search for a quasi-static balance occurring in the first instants, before 80 ms. The second indicator is the coupling behaviour while increasing (or decreasing) pressure loads, occurring whenever the acoustic pressure wave in the fluid meets the plate. This explains the choice of a water height H = 2m. The results are compared to the static theoretical solution proposed by [START_REF] Fourey | Développement d'une méthode de couplage fluide structure SPH Eléments Finis en vue de son application á l'hydrodynamique navale[END_REF]. The theoretical solution to this problem is not actually static since no viscosity and no damping in the material are present. The solution is expected to not diverge and to avoid any error amplification. The aim is thus to maintain the displacement of the midpoint of the plate within the envelop obtained for a constant loading case.

As a first step, a time step convergence study is presented for the CPS algorithm in Fig. 4.5. The first simulations are performed with the implicit Newmark scheme and without numerical dissipation (α s = 0). Note that Newmark and HHT schemes without dissipation lead to identical solutions. In all cases, the midpoint displacement is quickly damped at the early stage of the simulation. After t = 80 ms (first pressure wave return) the plate deformations become unstable. It can be observed that the oscillation amplitudes increase at each acoustic wave period (every 80 ms). Note that subdividing the time step is not sufficient to achieve a stable coupling in this case.

As shown in Fig. 4.6, adding a low numerical diffusion in the HHT scheme tends to cancel this noise and allows stabilizing the coupling. Similar vertical displacement time histories are obtained whatever the values adopted for α s , showing that this parameter does not introduce any phase shifting. After the transient solution, the α s parameter in the Newmark scheme increases exclusively the acoustic wave response peaks as α s increases. This result is counter-intuitive, and the origin of this effect should be further investigated. 4.7 shows the pressure field for different values of α s in the HHT scheme. It can be seen that the pressure field is stabilized whatever α s < 0. Note that the influence of α s is small for this case, pressure fields are similar in the range -0.05 > α s > -0.3.

Note also that the dissipation rate to introduce for stabilizing the coupling also depends on the material and the dynamics of the solid. Prior to launch the coupled simulation, a preliminary study on the structure vibration is needed to estimate the natural frequencies. This data is necessary to ensure that both solvers are designed with a sufficient dissipation and without altering the frequencies of interest. This is a key point for using this weak coupling strategy between SPH and FE methods.

All validations presented in FOUREY et al. ( 2017) have shown good agreement not only with the references, but also with the state of the art. Both CSS and CPS algorithms have been tested on various configurations. The CSS algorithm has revealed to be more stable than the CPS one, but in return CPS has shown better performances in terms of CPU time saving.

Discussion

A coupling strategy between SPH and FE methods has been presented in this chapter. A weak coupling has been chosen to couple these two methods for the sake of flexibility. The proposed coupling method can be used with any SPH or FE solvers and is relatively easy to implement. It appears that the present method is robust and compatible with an implicit time solution in the FE method, allowing larger time steps thanks to the absence of a prohibitive CFL condition. The next chapter is devoted to the analysis of the SPH-FE coupling from the energetic point of view, namely for a better understanding of numerical energy exchanges at the fluid/structure interface. Chapter 5

Energy considerations on the SPH-FE coupling strategy

In this chapter the energy conservation properties of our SPH-FE coupling strategy are investigated. According to [START_REF] Farhat | Two efficient staggered algorithms for the serial and parallel solution of three dimensional nonlinear transient aeroelastic problems[END_REF] and [START_REF] Piperno | Partitioned procedures for the transient solution of coupled aerobatic problems : Part II energy transfer analysis and three dimensional applications[END_REF] energy conservation and information transfers are primordial for coupling performances, accuracy and stability. In their work on coupling algorithms, MICHLER et al. ( 2004) underlined that partitioned algorithms do not permit the exact satisfaction of interface conditions, contrary to monolithic schemes. Any fluid-structure interaction model displays losses or gains of energy at the interface between fluid and deformable body. In their work, LI (2013), [START_REF] Li | A non-intrusive partitioned approach to couple smoothed particle hydrodynamics and finite element methods for transient fluid-structure interaction problems with large interface motion[END_REF] proposed a coupling strategy imposing an energy conservation at the interface. Conversely, FOUREY (2012) did not impose such a condition but succeeded in obtaining satisfactory results. [START_REF] Li | A non-intrusive partitioned approach to couple smoothed particle hydrodynamics and finite element methods for transient fluid-structure interaction problems with large interface motion[END_REF] did not make a coupling energy estimation analysis. Intuitively, one can think that in fast dynamics (as the aquaplaning phenomenon) a certain level of errors in the energy conservation would not be so important. Nevertheless, these errors could be significant for cases involving high pressure loading and velocities, even if the duration remains small. We consider that energy losses or gains should still be estimated in respect to other dissipative energies, to better understand the importance of enforcing conservation at the interface. Ensuring the rightful energy exchanges between the solid and the fluid is fundamental to build a correct model for fluid-structure interactions. Therefore the energy conservation properties of our SPH-FE coupling strategy are investigated. We first quantify the interface energy by following the same process as the one in [START_REF] Antuono | Energy balance in the δ-SPH scheme[END_REF]. Interface coupling energy is then compared to other energy components, especially with the numerical dissipation due to the stabilization of the SPH scheme and with the total fluid energy. Influence of various coupling parameters are also investigated. Finally various techniques for the determination of pressure loads on the deformable body are proposed and discussed. Here energy conservation properties are analysed by considering the case of a fluid interacting with a deformable body through SPH-FE coupling. The energy terms are expressed for two SPH schemes, with a particular focus on internal and boundary energies. The analysis start with the δ-SPH method, before making its analogy with the Riemann-SPH scheme actually used throughout this chapter. An energy balance over the entire domain is considered (fluid + solid), as the global domain must ensure the conservation of energy (Fig. 5.1). In the absence of outer contributions, several sources of energy creation or dissipation exist. The first source comes from the inside of each sub-domain. It corresponds to numerical diffusions and can be expressed as an internal irreversible energy. The second source of dissipation / creation comes from the treatment of the fluid structure interface. Determining the numerical diffusion in the SPH method is relatively straightforward, contrary to the FE method which requires to compare the energy at time t to its initial value. Note that the numerical diffusion in each sub-domain is negligible in comparison with other terms (at least for sufficiently converged solutions in these cases). In the following simulations, the numerical diffusion within the solid is considered as negligible compared to the fluid one. Indeed conservation properties are significantly better in the FE method than in the SPH one.

Energy considerations over the global domain

The second source of dissipation/creation is more subtle, as emphasized in this work. From the FE point of view the interface definition is quite obvious, as it corresponds to the wet surface of the deformable solid (Fig. 5.2) that is loaded with the fluid pressure. The energy E f luid/body transferred at the interface follows: However, to define the exact location of the interface is not obvious from the SPH method point of view. In a first approach, one can consider that it coincides with the wet surface (i.e. with the FE faces). Under this assumptions the energy transmitted at the interface could be the same from SPH (E body/f luid ) and FE (E f luid/body ) point of views. Nevertheless, the spatial operators of SPH make that the deformable body does not coincide exactly with the area of application of solid forces on the fluid. This area is located in the near boundary region (Fig. 5.2). In a second approach (more Hamiltonian) the interface could be considered as the whole set of particles contributing to solid boundary conditions (i.e. the particles having a support kernel truncated by a FE face). It is then possible to compute the power actually transmitted from the deformable solid to each fluid particle. The fluid-structure interface treatment is therefore quite subtle as the energy transmitted through the interface should be computed on different locations. Moreover this energy does not depend exclusively on pressure forces but also on other fluid characteristics as discussed further. At convergence the SPH and FE interface energies should tend to be the same:

E body/f luid ≈ E f luid/body , (5.2)

otherwise an interface coupling energy E Interf ace can be defined as:

E Interf ace = E body/f luid -E f luid/body .
(5.3)

ENERGY BALANCE IN THE δ-SPH SCHEME

The δ-SPH scheme assumes a constant mass for each particle i ( dm i dt = 0) with m i = ω i ρ i . In these conditions the variation of global kinetic energy within the whole domain writes: (5.12) where the artificial viscosity is considered as a pressure-like term

dE K dt = N F P i=1 ρ i ω i v i . g i - N F P i=1 N j=1 (P i + P j ) v i . ∇ i W ij ω i ω j + N F P i=1 N j=1 Π ij v i . ∇ i W ij ω i ω j ,
Π i = Π ij 2 .
Note that this equation can directly be obtained by multiplying the momentum equation by the velocity and then by summing over the fluid particles.

Internal energy

The variation of internal energy per unit mass E m I is expressed as:

ρ dE m I dt = -P div V + τ : D + r -div Q, (5.13)
where r is a heat source, Q a heat surface flux, τ : D the viscous dissipation, and -P div V a compressible power. In this study both r and Q are assumed to be zero, and an inviscid fluid is considered (τ : D = 0). Once more, the artificial viscosity is considered as a pressure-like term. This term has to be included in the internal energy equation since it acts as a pressure term for the variation of the particle volume. The variation of internal energy within the whole domain therefore writes:

dE I dt = - N F P i=1 N j=1 P i ( v j -v i ). ∇ i W ij ω i ω j + N F P i=1 N j=1 Π ij ( v j -v i ). ∇ i W ij ω i ω j . (5.14)
Let us focus on the first term of Eq. (5.14). Following [START_REF] Antuono | Energy balance in the δ-SPH scheme[END_REF] it is possible to distinguish two parts:

- 5.15) where P δ represents the power associated to the diffusive term in the mass conservation equation: (5.16)

N F P i=1 N j=1 P i ( v j -v i ). ∇ i W ij ω i ω j = N F P i=1 ω i ρ i P i dρ i dt dE C dt - N F P i=1 N j=1 P i ρ i ψ ij . ∇ i W ij ω i ω j P δ , (
P δ = N F P i=1 N j=1 P i ρ i ψ ij . ∇ i W ij ω i ω j ,
and E C is a compressible energy due to the compressible feature of the SPH method:

dE C dt = N F P i=1 ω i ρ i P i dρ i dt .
(5.17) CHAPTER 5. ENERGY CONSIDERATIONS ON THE SPH-FE COUPLING STRATEGY However, in a weakly-compressible SPH method this energy is expected to remain rather small.

It is now possible to express the internal energy into two parts E I = E C + E V as:

dE V dt = -P δ + N F P i=1 N j=1 Π ij 2 ( v j -v i ) ∇ i W ij ω i ω j .
(5.18) E V is an irreversible internal energy. It represents the energy lost to stabilize the SPH method. It does not affect the consistency of the SPH scheme since it tends to vanish as the spatial resolution increases.

Boundary energy

Considering solid boundaries, a certain amount of energy can be lost or gained by the fluid. This energy is called E body/f luid hereafter. Using the ghost particle method, [START_REF] Colagrossi | Numerical simulation of interfacial flows by Smoothed Particle Hydrodynamics[END_REF], the variation of the total energy verifies:

dE T OT dt = d(E P + E K + E C + E V ) dt = dE body/f luid dt , (5.19)
so that the summation of Eq. ( 5.9), Eq. ( 5.12) and Eq. ( 5.14) leads to: (5.20) where N g is the number of ghost particles inside the kernel support. Note that the wall velocity does not directly appear in this formulation: it is included within the ghost velocity v j . Using this approach, the transfer of energy from the solid to the fluid is estimated using all fluid particles having a truncated kernel near the boundary. Furthermore, nothing formally imposes dE body/f luid dt and dE f luid/body dt to be equal. There is indeed no direct correlation between Eq. ( 5.20) and Eq. (5.1). One can only expect that these two quantities are as close as possible. Furthermore considering fixed and flat boundaries without gravity leads to E body/f luid = 0. Indeed ghost particles have the same pressure as their mother particles, and opposite normal velocity. The conservation of the total energy is therefore ensured: dE T OT dt = 0.

dE body/f luid dt = - N F P i=1 Ng j=1 (P i + P j ) v i . ∇ i W ij ω i ω j + N F P i=1 Ng j=1 Π ij v i . ∇ i W ij ω i ω j - N F P i=1 Ng j=1 P i ( v j -v i ). ∇ i W ij ω i ω j + N F P i=1 Ng j=1 Π ij 2 ( v j -v i ). ∇ i W ij ω i ω j = - N F P i=1 Ng j=1 (P i v j + P j v i ). ∇ i W ij ω i ω j + N F P i=1 Ng j=1 Π ij 2 ( v j + v i ). ∇ i W ij ω i ω j ,

ENERGETIC STUDY OF THE SPH-FE COUPLING STRATEGY

Boundary energy

The boundary energy term is also similar:

dE body/f luid dt = - 1 γ i N F P i=1 Np j=1 (P i v j + P j v i ). ∇ i W ij S j ω i + 1 γ i N F P i=1 Np j=1 Π ij 2 ( v j + v i ). ∇ i W ij S j ω i .
(5.38)

Summary

This section was dedicated to energy considerations inside SPH schemes. It is now possible to compute the energy transferred at the solid interface (E body/f luid ). Consequently, E Interf ace can also be estimated for several test cases using Eq. ( 5.3). This study aims at analysing the amount of energy introduced or dissipated compared to the other energy components. Some differences with the energy balance proposed by [START_REF] Cercos-Pita | SPH energy conservation for fluid-solid interactions[END_REF] lead to the addition of an extra boundary term E ∂Ω V in Eq. ( 5.14):

dE ∂Ω V dt = N F P i=1 Ng j=1 Π ij 2 ( v j -v i ) ∇ i W ij ω i ω j .
(5.39)

In [START_REF] Cercos-Pita | SPH energy conservation for fluid-solid interactions[END_REF], the dissipated energy from the pressure-like term Π i = Π ij 2 does not involve the boundary conditions in the internal energy equation. However the model consistency is ensured: at convergence this term tends to disappear (Π ij → 0).

Energetic study of the SPH-FE coupling strategy

Two bi-dimensional test cases are considered here to analyse the energetic properties of the SPH-FE coupling. The hydrostatic water column on an elastic plate is not considered here due to the insignificant level of energy transferred at the interface compared to the others energy terms. Here the influence of the dissipation parameter α s inside the FE time integrator schemes (Newmark and HHT) is also investigated, as it rules the SPH-FE coupling stability [START_REF] Fourey | An efficient FSI coupling strategy between smoothed particle hydrodynamics and finite element methods[END_REF]. Note that only the CPS algorithm is used as it is the less stable according to Section 4.4. In a first approach, the ghost particle method is considered for the fluid boundary conditions inside the SPH solver.

Deformable beam impact

The first test case proposed is a violent free surface impact of a deformable beam structure driven at high velocity [START_REF] Scolan | Hydro-elastic behaviour of a conical shell impacting on a quiescentfree surface of an incompressible liquid[END_REF] similar to the one presented in [START_REF] Peseux | Hydrodynamic impact: Numerical and experimental investigations[END_REF], as illustrated in Fig. 5.3. Tab. 5.1 summarizes the material characteristics and the simulation parameters used. The beam is made of aluminium and impacts the free surface initially at rest with an inclination angle β = 10 degrees. The motion is imposed at the two extremities with STRATEGY C 0 1500.0 ms -1 CFL number 0.375 R/∆x SP H 3.0 Table 5.1: Physical and numerical parameters for the deformable beam impact. Fig. 5.4 provides the midpoint displacement time history of the beam for different spatial resolutions in the fluid. For ∆x SP H = 1 mm in the impact area, results are close to the semianalytical solution provided in [START_REF] Scolan | Hydro-elastic behaviour of a conical shell impacting on a quiescentfree surface of an incompressible liquid[END_REF], obtained by combining the hydrodynamic Wagner model and a linear model of elasticity for thin shells. Good agreements are also found with results from [START_REF] Li | A non-intrusive partitioned approach to couple smoothed particle hydrodynamics and finite element methods for transient fluid-structure interaction problems with large interface motion[END_REF] which enforced energy conservation at the fluidstructure interface. The beam deformation increases in the first instants, in accordance with the increase in pressure. It is nevertheless noted that the deformation of the plate in its midpoint finally appears slightly lower in the case of SPH-FE coupling than in the semi-analytical case. This difference can be explained by the non-consideration of the jet 5.4. ENERGETIC STUDY OF THE SPH-FE COUPLING STRATEGY in this solution. It can also be explained by an underestimation of the fluid pressure in the near body area according to Fig. 5.6. Nevertheless, note that the peak occurrence is correctly captured for ∆x SP H = 1 mm, as well as the pressure decrease after the peak. Moreover, a good agreement is observed for the vertical force applied onto the deformable structure for ∆x SP H = 1 mm (see Fig.

CHAPTER 5. ENERGY CONSIDERATIONS ON THE SPH-FE COUPLING STRATEGY

From the energetic point of view, an energy transfer occurs from the deformable beam to the water (Fig. 5.7). This energy results in fluid particle motions, visible through an increase of mechanical energy E M (E M = E P + E K ). Water particles are propelled out of the beam extremity as a high-speed free surface jet (visible in Fig. 5.3). E C converges fast (see Fig. 5.8). The compressible nature of the SPH scheme is used to correctly model this impact case. The fluid is first highly pressurized and then depressurized in the jet. The maximal velocity reaches 300m.s -1 , corresponding to a Mach number of 0.2. C 0 has therefore to be chosen as equal to the actual water sound speed. This is necessary to correctly estimate the part of the energy transferred from the deformable beam which is converted into compressible energy ( E C dt > 0) instead of mechanical energy, before being released at the end of the simulation ( E C dt < 0). The compressible effects are not of major importance but cannot be neglected, so that an incompressible method would miss a part of the flow physics. Note that a significant amount of energy is also dissipated to maintain the stability in the SPH fluid domain. According to the consistency property it decreases as the spatial resolution increases. However it tends to mask the capture of compressible effects, even with the finest spatial resolution. Nevertheless, in this configuration E C + E V is about 10% of E K . The energy errors from E V remain therefore acceptable using ∆x SP H = 1mm. Concerning the coupling interface energy, Fig. 5.9 shows that the fluid receives less energy with the ghost particle method (E body/f luid ) than the energy lost by the aluminium beam through pressure forces (E f luid/body ). Following Eq. ( 5.3), E Interf ace is negative. Note that dE body/f luid dt is strictly equal to dE T OT dt according to Eq. (5.19). Fig. 5.10 shows the energy variation in the global domain ∆E F +S = (E T OT + E S ) -(E 0 T OT + E 0 S ) compared to the energy yielded from the outside E Outside , where E S is the solid total energy deduced 5.4. ENERGETIC STUDY OF THE SPH-FE COUPLING STRATEGY from Code Aster outputs. E Outside corresponds to the energy brought to impose a constant velocity (30 m.s -1 ) on the beam, which is estimated using the FE solver outputs. The global system loses energy through the fluid-structure interface. This phenomenon can lead to severe problems. Indeed, the accumulation of energy errors can modify the accuracy and the stability of the solution. Note that Fig. 5.9 and Fig. 5.10 are very similar. This was expected since E Outside mostly depends on the energy recovered by the water from the FE point of view (i.e. E f luid/body ) in order to maintain a constant velocity (30 m.s -1 ) on the beam. Furthermore, ∆E F +S represents the energy variation in the global domain, which mainly corresponds to the fluid energy variation (i.e. E body/f luid ). Only a small part is taken over by the beam principally through elastic energy. Note that E f luid/body displays a lower convergence rate than E body/f luid , so that E f luid/body requires a higher spatial resolution to converge. However the adopted spatial resolution is sufficiently high to correctly simulate this deformable beam impact (Fig. 5.5). Furthermore, the consistency of our model seems ensured as both E f luid/body and E body/f luid tend to converge towards the same value. Strong variations of E Interf ace occur during the impact, as visible in Fig. 5.11. From t = 0 ms to t = 1.5 ms, E Interf ace represents more than 20% of E T OT . At the beginning of the impact errors are close to be equal to E T OT for ∆x SP H = 1 mm. For coarser fluid spatial resolutions, the maximal ratio is lower because E T OT is higher. However, Fig. 5.9 shows that the total energy is very small at this stage (order of 10 4 J) compared to the following (order of 10 5 J). The simulation remains therefore acceptable despite of these important errors, at least for the highest spatial resolution. ∆x SP H = 1 mm is clearly just sufficient, coarser resolutions are not acceptable. For the finest resolution, the energy ratio quickly decreases under 0.1 after t = 1.5 ms. It finally converges towards a value close to 7% of E T OT . At this time E Interf ace has the same energy level as E V (see Fig. 5.12). The lack of accuracy in the energy conservation over the global domain is not acceptable CHAPTER 5. ENERGY CONSIDERATIONS ON THE SPH-FE COUPLING STRATEGY regarding the other terms. Note also that Fig. 5.12 shows that varying the parameter α s in the FE integration scheme has no influence on the solution for α s < -0.05 with the Newmark scheme. However in practice α s ≤ -0.1 is needed for stability purposes using the HHT scheme. Instabilities appear for lower values of α s [START_REF] Fourey | An efficient FSI coupling strategy between smoothed particle hydrodynamics and finite element methods[END_REF] leading to divergences of E Interf ace . Fig. 5.14 confirms the above statements, by presenting snapshots of the computed solution at various instants with different FE time integrator schemes. These snapshots are plotted with pressure and Von Mises stress contours respectively for the fluid and the structure. Strong non physical perturbations are observed in the pressure field for HHT with α s = -0.05 (especially after t=2ms), leading to strong vertical oscillations as visible in Fig. 5.13. In this configuration, increasing the stabilization parameter α s helps in maintaining a stable coupling and reducing the energy errors (Fig. 5.12).

Dam-break flow through an elastic gate

A case involving relatively low dynamics is considered here. The fluid initially at rest is contained within a rigid tank closed using an elastic gate on the left side. At t = 0 s the elastic gate is released to let the water escape. This test case has been first introduced by ANTOCI ( 2006), according to the configuration presented in Fig. 5.15 and Tab. 5.2. In their study, ANTOCI et al. (2007) used a monolithic approach with a SPH-SPH coupling to model this FSI problem and proposed some comparisons with their experimental data. In the present work, fluid particles are distributed uniformly, and finite elements are used to model the elastic gate. The FE mesh is composed of 4 width-wise and 40 lengthwise elements. The behaviour of rubber is considered as incompressible with a Poisson coefficient close to 0.5, and the non-linear strain-stress curve obtained experimentally is used to characterize its behaviour (Fig. Fig. 5.17 compares the experimental and numerical time histories of the gate tip displacements for a stabilization coefficient α s = -0.3 in the Newmark scheme. Numerical results are in good agreement with the experiments from [START_REF] Antoci | Numérical simulation of fluid structure interaction by SPH[END_REF], even with the coarsest resolution (∆x SP H = 1 mm). In particular, the global trend of the gate deformation time history is correctly captured. A small offset can be noticed but it tends to decrease as the spatial resolution increases. The abrupt decrease of the displacement observed experimentally from instant t = 0.32s is not reproduced using our coupling, but it can be attributed to a leakage of fluid between the tank and the gate visible on the experimental snapshots. Numerical results are also in good agreement with [START_REF] Li | A non-intrusive partitioned approach to couple smoothed particle hydrodynamics and finite element methods for transient fluid-structure interaction problems with large interface motion[END_REF], where the energy conservation is enforced at the interface. The differences observed can be explained by small discordances in the rubber material law. From the energetic point of view, the water transfers its energy to the gate which is converted into elastic energy (Fig. 5.18). The flow is purely incompressible (E C is negligible). In the simulation, the initial water potential energy is thus mainly converted into kinetic energy (Fig. 5.19), but a small part is lost through the boundary conditions (Fig. 5.20). Note that E S , E P and E K are converged already with the coarsest fluid spatial resolution. The energy lost through fluid boundary conditions (exactly equal to the variations of E T OT ) is lower than the energy provided to the elastic gate by pressure forces (see Fig. 5.20). Consequently, an unexpected additional energy is introduced at the body interface. E Interf ace is therefore positive and the total energy is slightly increased. Fig. 5.21 shows the energy variation of the global domain

E F +S -E 0 F +S E 0 F +S
which increases while no energy is introduced from the outside. However the consistency of our model seems ensured since

E F +S -E 0 F +S E 0 F +S
tends towards zero as the spatial resolution increases. 2017), a small dissipation modifies the coupling stability, so that some instabilities tend to appear at low dissipation levels in the solid. For HHT coefficient α s = -0.1, E Interf ace increases much faster than for higher stabilization levels. Fig. 5.23 shows the displacement of the gate tip for various values of the stabilization parameter and confirms that the displacement is strongly erroneous with the HHT scheme using α s = -0.1. Instabilities in the structure solution lead to high frequency displacements of the interface, and are responsible for non physical acoustic waves within the fluid (Fig. 5.24).

All other cases plotted in Fig. 5.22 show negligible energy errors (less than 1%) when compared to E T OT . Nevertheless, as shown in Fig. 5.25 E Interf ace still displays strong errors with respect to E V even using the finest spatial resolution. At the beginning of the simulation, errors are close to be equal to 18 times E V (see Fig. 5.26). However, E V is very small at this stage (order of 10 -3 J) compared to the following (order of 10 -2 J). Then this ratio quickly decreases under the value 2 for t > 0.15s. As in Fig. 5.12, E Interf ace and E V have the same level. An interpretation is that the numerical dissipation needed to stabilize the SPH computation mainly occurs near to the deformable body where the flow is the most perturbed. As a result, E Interf ace may be mostly rules by the artificial viscosity term. However this assumption would need a further investigation. These simulations therefore remain acceptable although large apparent errors appear at the first instants. Nevertheless these errors are unacceptable in terms of accuracy about the energy conservation. Furthermore, E Interf ace tends to diverge for low numerical diffusion levels while using the HHT scheme, according to the stability criterion describes in [START_REF] Fourey | An efficient FSI coupling strategy between smoothed particle hydrodynamics and finite element methods[END_REF]. STRATEGY

Coupling improvements

Fluid-structure loading

The previous section emphasized the non conservation of interface energy and outlined the need for improving the coupling scheme proposed. A correction criterion or a new loading method is needed to overcome these limitations and to improve the coupling robustness and accuracy. For algorithmic simplicity and genericity purposes, we chose to avoid any hard-coding within the FE solver. The improvements proposed address exclusively fluid pressure loads provided to the FE solver and interface energy conservation.

Force conservative formulation

Our first approach consisted in determining fluid pressure loads as an average of neighbourhood pressure performed through Eq. (4.1). The second approach proposed here resides in ensuring reciprocal forces between fluid and solid media. This is achieved through the summation of ghost particle interactions with fluid particles, as proposed for instance in [START_REF] Doring | Développement d'une méthode SPH pour les applications à surface libre en hydrodynamique[END_REF][START_REF] Bouscasse | Nonlinear water wave interaction with floating bodies in SPH[END_REF] for rigid bodies. The force can be expressed as:

F f luid/body = -F body/f luid = - N F P i=1 Ng j=1 [-P i -P j + Π ij ] ∇ i W ij ω i ω j , (5.40)
Obviously in the case of a deformable body this force should be decomposed on each panel k, so that:

F f luid/body = Np k=1 F k f luid/body , (5.41)
with

F k f luid/body = - N k i=1 Ng j=1 [-P i -P j + Π ij ] ∇ i W ij ω i ω j , (5.42)
where N P is the number of panels and N k is the set of particles seen by the panel k (see Fig. 4.2). As the calculation of ghost interactions is needed anyway for imposing boundary conditions, pressure forces F k f luid/body are transferred to the FE solver at the end of each first Runge-Kutta stage. Using NFM, the principle is the same but using the summation of panel interactions with fluid particles.

Energy conservative formulation

An energy conservative exchange at the interface is also considered. At each instant, it is possible to compute the boundary energy E body/f luid of the previous time step. A corrective coefficient can be introduced in the force F f luid/body calculation as:

F f luid/body = -(1 + ) N F P i=1 Ng j=1 [-P i -P j + Π ij ] ∇ i W ij ω i ω j , (5.43)
The new variations of the boundary energy transmitted from the fluid to the body STRATEGY E f luid/body can be expressed as:

dE f luid/body dt = N P k=1 -(1 + ) N k i=1 Ng j=1 [-P i -P j + Π ij ] ∇ i W ij ω i ω j . v k = (1 + ) dE f luid/body dt ,
(5.44) On the contrary the variation of E body/f luid remains unchanged since it does not depend on . The coefficient is chosen to impose the same evolution for E f luid/body and E body/f luid leading to: According to the time step, E f luid/body and E f luid/body can be very different from each other. An average procedure is thus performed on the previous time steps to avoid force discontinuities and a limitation value is imposed so that:

E f luid/
-0.15 < < 0.15 (5.47) This can be justified by the fact that the correction remains small provided that the energy conservation is maintained.

Numerical validations

All the following results are performed using the ghost particle method, but similar tendencies are observed using the NFM method (see Appendix A).

5.5. COUPLING IMPROVEMENTS

Deformable beam impact

The deformable beam impact is performed here with the above improvements and using the same spatial fluid resolution ∆x SP H = 1mm. In the last section coupling errors tended to be higher than 7% of E T OT using the averaging pressure technique (see Fig. 5.11). Using the force conservative formulation, E Interf ace is still negative (the global system still lose some energy through the fluid structure interface), but energy errors converge towards a value close to 2.8% of E T OT (see Fig. 5.27). The accumulation of energy errors is therefore reduced thanks to the force conservation formulation. Another advantage of this approach is to avoid any overestimation of forces in cases where only few particles are located in the near panel area. Only ghost particle forces on the panel k are considered, and pressure loads are therefore more accurate. As a consequence the ratio between E Interf ace and E V is also strongly decreased (Fig. 5.27), as it is now close to 0.4 instead of 1.4 in Fig. 5.12.

Figure 5.27: Time history of ratios between E Interf ace and E T OT (left) and between E Interf ace and E V (right) for different stabilization coefficients, force conservative formulation.

Fig. 5.28 presents the time history of the same energy ratios, but obtained this time with the energy conservative formulation. Energy ratios converge towards zero at the end of the simulation and E Interf ace is not strictly equal to zero due to the limitation imposed on the corrective coefficient . Nevertheless it quickly converges towards zero in all cases, emphasizing the benefits of such improvements. Note also that the E Interf ace > 0 here, the corrective coefficient seems underestimated. Nevertheless the origin of this effect should be further investigated.

The robustness property of the SPH-FE coupling strategy is therefore strongly improved. Indeed the SPH-FE coupling is more stable for lower dissipative coefficients α s in the HHT scheme: even HHT with α s = -0.05 provides correct results, contrary to the first implementation. E Interf ace is no more affected by the dissipation in the structure part. Varying the parameter α s does not have a significant influence on the solution behaviour. 5.5. COUPLING IMPROVEMENTS Fig. 5.30 presents some snapshots of the computed solutions for the various formulations discussed in this paper. These snapshots are plotted with pressure and Von Mises stress contours respectively for the fluid and the solid. Using the averaging pressure procedure, the numerical diffusion introduced by the HHT is not sufficient. As a result the high frequency displacement of the deformable body interface generates some acoustic waves in the fluid. Using the force and energy conservation procedures proposed here, pressure loads are more accurate and regular. Less diffusion is therefore needed to obtain a robust coupling. The pressure field is slightly more regular with the energy conservation procedure, especially at t = 2ms as shown in Fig. 5.30. The results using NFM are available in Appendix A, which are very similar to these ones. All previous results show some small differences between the two boundary conditions.

Dam-break flow through an elastic gate

The dam-break flow through an elastic gate is performed using a spatial resolution ∆ SP H = 0.25mm in the fluid. Here again, the accumulation of energy errors has been strongly reduced using the conservation improvements proposed (Fig. 5.31 and Fig. 5.32). It is particularly visible on the ratio between E Interf ace and E V .

Both formulations increase the coupling robustness regarding the lowest dissipation levels in the FE time integrator scheme. Numerical stability is maintained even at low stabilization coefficients (Fig. 5.33). Force and energy conservative formulations show good agreements with the experimental data, but also with the numerical results from [START_REF] Li | A non-intrusive partitioned approach to couple smoothed particle hydrodynamics and finite element methods for transient fluid-structure interaction problems with large interface motion[END_REF] which enforce energy conservation at the fluid-structure interface using a very different approach. This test case shows two main phases: from t = 0 s to t ≈ 0.14 s the fluid yields its energy to the elastic gate, Fig. 5.34 (i.e the elastic gate is pushed). Then the solid gate gives back a part of its energy to the fluid from t ≈ 0.14 s to t = 0.4 s. Regarding the force time history, both formulations provide slightly different results (see Fig. 5.35). During the first phase, the dynamics of the flow are the same. E Interf ace is very small and both formulations provide similar results. In the second phase, some differences appear as the gate comes back. Using the energy conservative formulation, the elastic gate tends to come back more slowly (Fig. 5.34), as the force is increased to be in energy agreement. Fig. 5.36 shows some snapshots of the computed solution for the different formulations proposed. It highlights the benefits of the new loading formulations compare to the averaging pressure procedure. High frequency pressure waves in the fluid domain are avoided thanks to the gain in accuracy on the pressure loads. Here again the last two formulations provide very similar results.

Force conservative formulation

Energy conservative formulation

Discussion

The analysis of a SPH-FE coupling method has been provided from an energetic point of view with emphasis on interface coupling energy. An energy balance study has been proposed for both δ-SPH and Riemann-SPH schemes, highlighting the presence of a term representing the amount of energy dissipated or created at the fluid-structure interface. Inaccurate fluid loading leads to an accumulation of energy errors which are responsible for coupling instabilities. In this chapter two different fluid loading procedures have been proposed to increase the conservation property of the coupling. The energy errors have been decreased and the robustness properties of the coupling have been improved. Only small differences were observed between the force and energy conservative formulations proposed. However the energy conservative procedure requires the energy term calculation for all particles, increasing the computational costs. As a consequence, the force conservative formulation can be preferred regarding the compromise between coupling accuracy, robustness and CPU time. Furthermore, two different fluid boundary conditions techniques (ghost particle method and NFM) have been considered, showing equivalent results in each case.

Chapter 6

Computational improvements of the SPH-FE coupling strategy

The energy conservation properties of our SPH-FE coupling have been investigated in the last chapter. This study allowed us to validate the FSI strategy from an energetic point of view, but also to improve its robustness and accuracy properties. Nevertheless, the features related to the computational cost were not analysed. This thesis aims at studying the aquaplaning phenomenon in an industrial context. Thus, the SPH-FE coupling strategy should be optimized to deal with 3D complex flows. Depending on the physics of the FSI problem, one of the two sub-domains needs smaller time steps. For the present SPH-FE coupling, time steps are expected to be smaller within the fluid domain: the CFL condition imposes small fluid time steps due to the fine spatial resolution. Furthermore, computational times are larger in the serial FE solver given the parallel nature of SPH-flow OGER et al. (2016a). Contact interactions between the tire and the road also require fine solid spatial resolutions, inducing large CPU times for the FE simulation. Consequently, CPU times are mostly ruled by the solid part. This chapter is dedicated to proposing and analysing some improvements to reduce the computational costs.

Two complementary approaches are considered. The first approach aims at reducing the FE solver involvement in the coupling computational time. This approach is based on a new coupling algorithm, permitting a reduced number of calls to the FE solver. This procedure is possible thanks to the large ratio between solid and fluid time steps. The second approach proposes to use the local particle refinement defined in Section 2.8.2 to decrease the CPU cost on the SPH side. A linear extrapolation of the structure displacement is carried out for each fluid step to according to:

x t n+l i = x t n i + l∆t v t n i , ∀l = 1, ..., p (6.1)
where x i and v i are the positions and velocities of the node i. This is necessary to avoid fluid pressure wave generation at the fluid-structure interface.

In addition, a regulation of the ratio p is performed according to the variations of the fluid loads at the end of each loop:

           p = p -1 if 0.075 < k F k,n+p f luid/body -k F k,n f luid/body k F k,n
f luid/body ≤ 0.15

p = p + 1 if k F k,n+p f luid/body -k F k,n f luid/body k F k,n f luid/body ≤ 0.075 p = 1 otherwise , (6.2)
When the force variations reach 7.5% between two loops, the assumption F n+p f luid/body ≈ F n f luid/body (described in Section 4.3) is not accurate. This procedure allows us to increase the number of fluid-structure interactions when it is needed, to correctly capture all coupled effects. Furthermore if a variation exceeds 15%, the ratio is directly reduced to 1. Otherwise the ratio p is increased by one, within a limit defined at the beginning of the simulation.

Furthermore, each loop should ensure:

∆t F E = p∆t SP H , (6.3)

With the CPS procedure, it is theoretically impossible to define the SPH time step size beforehand (∆t n+1 , ∆t n+2 ,...) due to the variations of the sound speed. Therefore, a prior estimation of the time step is used, by choosing the time step at instant t n imposed by the CFL condition and by limiting its size for safety purpose as: ∆t n+i = 0.8∆t n , ∀i = 1, ..., p (6.4)

∆t F E = 0.8p∆t n . (6. (6.9)

Numerical validations

This method has been tested on two different tests cases. The aim is to validate this approach and to estimate the maximum p ratio maintaining satisfactory results. All the following simulations are performed using the Newmark scheme with α s = -0.3 for the solid time integration, as it is the most stable. Several values of the time step ratio are studied here. The ghost particle method is used for the solid boundary conditions in the SPH solver. The results with the NFM give exactly the same trends, for that reason they are not presented here. Numerical results are compared regarding fluid forces, body displacements, fluid pressure fields, etc.

Deformable beam impact

The deformable beam impact is performed using a fluid spatial resolution ∆x SP H = 1 mm and with the simulation parameters described in Tab. 5.4.1 and Fig. 5.3. Time step ratios p from 1 to 20 have been tested. Fig. 6.2 and Fig. 6.3 show the time histories of the vertical fluid force and of the beam midpoint relative displacement. Very good agreements are observed, even for large ratio values (i.e. p = 20). It was expected since the time steps are very small in the fluid domain. Indeed, the numerical sound speed has been set to 1500 m.s -1 in the water (Tab. 5.1) to consider the physical fluid compressibility. In contrast, the solid has its first three modes of vibration at 575, 1526and 2864Hz according to FOUREY (2012). The vibration are therefore correctly captured even using p = 20. Fig. 6. 4 and Fig. 6.5 show some snapshots of the computed solution for different time step ratios using the CPS procedure. These snapshots are plotted with pressure and Von Mises stress contours for the fluid and the solid respectively. Identical trends are obtained with the CSS procedure (see Appendix B.1). The pressure field is preserved when the time step ratio is increased from 1 to 16 and no irregularity or instability is observed. CHAPTER 6. COMPUTATIONAL IMPROVEMENTS OF THE SPH-FE COUPLING STRATEGY While using a ratio p = 20, some small perturbations begin to appear in the pressure field (see Fig. 6.6). Some streaks can be observed in the flow. These non-physical pressure waves come from the actualisation of the node positions. Indeed, using large p ratio values the kinematic extrapolations Eq. (6.1) generate some significant errors. The differences between the last extrapolated positions at t n+p-1 and the newly computed position at t n+p are too large. As a result, a pressure wave is created at the fluid-structure interface. The perturbation then propagates in the fluid domain, leading to a loss of accuracy and stability. As a consequence, the p ratio cannot be increased too much to maintain the stability and accuracy properties. This observation underlines the importance of the extrapolation Eq. (6.1). The absence of such an extrapolation leads to strong errors, even at lower p ratio values as shown in Fig. 6.7. In that case, the vertical fluid force and the beam relative displacement are overestimated (Fig. 6.8). Note that linear extrapolation seems sufficient, as second order extrapolation has also been tested without any significant improvement, even for p = 20. CHAPTER 6. COMPUTATIONAL IMPROVEMENTS OF THE SPH-FE COUPLING STRATEGY

Dam-break flow through an elastic gate

The dam-break flow test case is studied using the fluid spatial resolution ∆ SP H = 0.25 mm and with the same parameters as described in Section 5.4.2. This test case is naturally less stable than the previous one. Acoustic pressure waves generated at the deformable body interface are reflected onto the opposite wall. These reflections lead to high frequency fluid loading, which tends to destabilize the fluid computation as described in [START_REF] Fourey | An efficient FSI coupling strategy between smoothed particle hydrodynamics and finite element methods[END_REF]. The dissipation inside the solid time integrator scheme maintains the coupling stability. However the pressure waves generated by the differences between the newly computed positions at t n+p and the last extrapolated positions at t n+p-1 affect the solution. The number of back-and-forth travels of the pressure wave is directly linked to the sound speed. As a result, C 0 has been reduced to 25 m.s -1 to minimize the number of wave reflections, so as to reinforce the numerical stability.

Fig 6.9 shows the time history of the horizontal force. A CPS procedure with p > 8 generates large high frequency oscillations, which degrade the fluid loading quality. The coupling is clearly not stable, as visible in the pressure field (see Fig. 6.10). Nevertheless, very good results are obtained for ratios p≤4 with both algorithms. Note that the ratio p = 16 is in agreement with the natural vibration of the gate, the first mode is at 11 Hz according to [START_REF] Fourey | Développement d'une méthode de couplage fluide structure SPH Eléments Finis en vue de son application á l'hydrodynamique navale[END_REF]. Theses force oscillations are not observed using the CSS procedure, which is more stable due to the consistent imposition of loads in time. As it could be expected, the maximum value allowed for p depends on the stability property of the coupling. Fig. 6.12 and Fig. 6.13 provide some snapshots of the computed solution for different time step ratios using the CSS procedure (the rest of the results with the CPS procedure are available in Appendix B.2). These snapshots are plotted with pressure and Von Mises stress contours for the fluid and the solid respectively. Good agreements are observed between ratios from 1 to 8. Larger differences begin to appear with the ratio p = 16. Acoustic effects are more visible, but are still acceptable. The velocity field remains unchanged.

Summary

The possibility to reduce the FE solver involvement in the total computational cost has been proved in this section. Two different coupling algorithms have been proposed to minimize the number of FE time steps by imposing different time steps to each sub-domain. It prevents the time step requirements of the solid sub-domain from being imposed by the fluid sub-domain. This procedure is possible since the explicit fluid computation requires smaller (in general case) time steps than the solid one. For each case, the maximum ratio p is imposed as a compromise between stability features and the computational times. The CSS procedure outperforms the CPS one this regard, due its consistent imposition of the fluid loading in time. 6.2. ADAPTIVE PARTICLE REFINEMENT (APR)

Adaptive particle refinement (APR)

The adaptive particle refinement technique is an additional solution to reduce the CPU times. Here, we proposed to test the adaptive particle refinement proposed by [START_REF] Chiron | Analysis and improvements of adaptive particle refinement (APR) through CPU time, accuracy and robustness considerations[END_REF] (described in Section 2.8.2) in the SPH-FE coupling context. Several APR configurations are studied. Numerical results are then validated regarding fluid forces and body deformations, etc.

Deformable beam impact

The deformable beam impact is simulated here using the particle refinement instead of the variable spatial resolution. Several refinement boxes are defined according to Fig. 6.14. The last box corresponds to the finest resolution zone (∆x SP H = 1 mm). Here, six refinement levels are used leading to ∆x SP H = 3.2 cm in the coarsest resolution zone. The number of particles is divided by five compared to a fully refined simulation at t = 0 s. In the first attempt, all refinement boxes are fixed. Very good results are observed concerning the fluid vertical force (Fig. 6.15) and beam midpoint relative displacement (Fig. 6.16). Coupling accuracy and stability are not influenced by the number of refinement levels whatever the number of refinement boxes used.

However, fixed boxes restrain the remarkable possibilities of the adaptive particle refinement. Indeed in the application targeted in this work, the tire will cross a puddle of water in an aquaplaning simulation, so that the zones of interest will not be fixed. Upstream particles will not need to be as fine as the ones near the tire. Therefore, in a second approach we made refinement boxes follow the body movements at each time step. Good agreements are observed when compared to fixed boxes (Fig. 6.17 and Fig. 6.18). Nevertheless, note that small deviations are observed for both fluid force and body displacement. An explanation STRATEGY In these configurations the CPU time has been reduced by 50 %. The efficiency of the APR technique is directly linked to the number of particles by processor, [START_REF] Chiron | Couplages et améliorations de la méthode SPH pour traiter des écoulements multi-échelles temporelles et spatiales[END_REF]. A low ratio between the number of fluid particles per processor has been considered here, explaining the low gain in the computational costs.

Dam-break flow through an elastic gate

The dam-break flow through an elastic gate is also simulated using the particle refinement technique to limit the number of involved particles. Several boxes are defined according to Fig. 6.21. The particle spacing required between two boxes prevents us from using more than four refinement levels due to the buffer zones. The fluid resolution corresponds to ∆x SP H = 0.25 mm in the finest zone. Using four levels, the resolution reaches ∆x SP H = 2 mm in the coarsest zone. The number of particles is divided by six compared to the fully refined resolution at instant t = 0 s. All refinement boxes are fixed. In this configuration the CPU time has been reduced by 40 %. Once more, the gain seems limited by the low number of fluid particles per processor. Here again, very good results are observed concerning the fluid horizontal force (Fig. 6.23) and the gate tip displacement (Fig. 6.22). No instability appeared in the SPH-FE coupling calculations.

APR boxes have another valuable advantage for SPH-FE coupling models. They behave as low-pass filters onto the SPH fluid loading. High-frequency pressure waves generated at the interface are partially dissipated through each APR box crossing. The coarsest resolutions are not sufficient to capture these highest frequencies. Therefore, this low-pass filter naturally increases the SPH-FE coupling stability. It is then possible to look for higher p ratios. Some simulations can be performed using the CPS algorithm with p = 8 (Fig. 6.24) contrary to the previous section (recalled in Fig. 6.25).

This behaviour can also be observed on the fluid pressure field. Fig. 6.26 provides some snapshots of the computed solution for various APR configurations. These snapshots are plotted with pressure and Von Mises stress contours for the fluid and the solid respectively. Good agreements are observed between the different each configurations. CHAPTER 6. COMPUTATIONAL IMPROVEMENTS OF THE SPH-FE COUPLING STRATEGY

Discussion

This chapter focused on improvements of the coupling regarding the computational costs. Two different approaches, which can be combined, have been proposed to reduce the computational costs of our simulations. The use of different time steps in the fluid and solid sub-domains provides a good compromise between CPU time, accuracy and stability. On the other side, it has been proved that the APR allows reducing the SPH CPU times by minimizing the number of considered particles without decreasing the solution accuracy. On the contrary, the APR acts as a low-pass filter of the fluid loading providing some benefits. The refinement method can therefore increase the coupling stability in some cases. These two approaches can be used simultaneously to decrease the computational time for a given CPU resource. A good compromise among stability, accuracy and computational times has therefore been obtained using these new tools. ∆x F E et al. (2007) and the 2D results of [START_REF] Li | A non-intrusive partitioned approach to couple smoothed particle hydrodynamics and finite element methods for transient fluid-structure interaction problems with large interface motion[END_REF]. The global trend of the gate deformation is captured. The abrupt decrease of the displacement observed experimentally from instant t = 0.32 s is still not reproduced by our coupling. No leakage of fluid is observable in the simulation contrary to the experiment. This was expected since free-slip condition is considered between the gate and the tank. Furthermore, the 3D numerical results are in good agreement with the 2D results which as expected since this problem with free-slip condition on the lateral wall is quasi-2D. Fig. 7.3 compares the experimental and numerical water height time histories just behind the gate (x = 0.4 m) and in the middle of the reservoir (x = 0.45) m. Good agreements are observed despite a small underestimation of the water level. Exactly as for the 2D results, it seems that the present results show a small persisting offset of about of 0.04 s. The origin of this difference is not clearly identified, but may be imputed to experimental uncertainties in the triggering, since the experimental water level tends to remain stable in the first instants while the gate has started to move. Fig. 7.4 shows the fluid pressure field at various instants for the 2D and 3D models and The following aquaplaning simulations will be performed using Abaqus software as FE solver. Nevertheless all previous simulations were performed with Code Aster. Consequently, a comparison is made between these two FE solvers to check the robustness of the SPH-FE coupling strategy. The APR method is also used to reduce the number of involved particles, as 3D simulations are demanding in terms of CPU resources. The following simulations are performed with a single APR box according to the configuration presented in Fig. 7. 5. With one refinement level, the fluid spatial resolution reaches ∆x SP H = 2 mm in the coarser zone, while the finer resolution is still maintained as ∆x SP H = 1 mm.

Sloshing tank interactions with an elastic beam

The second and third validation cases are performed on a sloshing tank interacting with a deformable clamped beam. These test cases have been introduced by IDELSOHN et al. ( 2008) and described in BOTIA-VERA (2015). The initial configurations are provided in Fig. 7.9. A rectangular tank containing a clamped beam is filled with sunflower oil (two different levels). A forced roll motion is then applied to the tank, as shown in Fig. 7.10. The center of rotation is located at the midpoint of the tank bottom. Furthermore, the period is chosen so as to coincide with the critical sloshing frequencies, according to the filling level. Initially another test case concerning a hanging beam has also been proposed.

It has not been investigated here due to a lack of time.

Mid depth oil

Shallow depth oil Two blowholes are located at the top of the tank to let the air circulate freely. The air phase can therefore be legitimately omitted provided that no air entrapment occurs during the evolution. The particle refinement is used to limit the number of fluid particles. The viscous nature of the flow is modelled using the Monaghan and Gingold formulation expressed in Section 2.6.2. The FE mesh is composed of 8 width-wise, 8 depth-wise and 80 length-wise elements. The beam material is a commercial neoprene rubber. The experimental strain-stress curve is not available here. The rubber behaviour is considered as incompressible with a Poisson coefficient close to 0.5 and the beam deformation is modelled by a linear elastic law, with E = 6 MPa. In each case, a gap of 2.9 mm exists between the clamped beam and the tank in the transverse z direction (see Fig. 7.11). The impact of these gaps has not been investigated. The fluid spatial resolution used is not sufficient to correctly capture the flow in these areas, but it could be the topic of a further work. Finally, the boundary conditions are modelled using ghost particles.

Mid depth oil

Shallow depth oil

Clamped beam in mid depth oil configuration

In this case, the tank is filled with 114.8 mm oil depth. This test case has already been studied in 2D by [START_REF] Yang | Free-surface flow interactions with deformable structures using an SPH-FEM model[END_REF] and [START_REF] Liao | A coupled FDM-FEM method for free surface flow interaction with thin elastic plate[END_REF]. Nevertheless, the flow concerns a viscous fluid in a narrow tank. The effects of the no-slip condition on the lateral walls cannot be neglected. Consequently, the use of a 3D model is crucial here. Three refinement levels are considered, according to the configuration in Fig. 7.12, where ∆x SP H = 4 mm in the coarsest zone. All simulation parameters are presented in Tab. 7.2.

Clamped beam in shallow oil configuration

The tank is now filled with 57.4 mm oil depth. A 3D model is still considered. 2D numerical results can be found in [START_REF] Degroote | Partitioned simulation of the interaction between an elastic structure and free surface flow[END_REF][START_REF] Paik | Fluidstructure interaction for an elastic structure interacting with free surface in a rolling tank[END_REF][START_REF] Paik | Fluidstructure interaction for an elastic structure interacting with free surface in a rolling tank[END_REF]. Three refinement levels are considered here according to the configuration displayed in Fig. 7.17, ∆x SP H = 4 mm in the coarsest zone. All simulation parameters are presented in Tab. 7.3.

The amplitude of the forced roll motion is also about 4 degrees (Fig. 7.18). This test case is more complicated than the previous one since the beam deformations are relatively small. In this configuration, the maximum beam deformation does not exceed 10 mm (see Fig. 7.19). The obtained numerical deformation is comparable to the experimental one, both in phase and in amplitude. Nevertheless, the numerical results are farther from the experimental solution than those for the mid depth case (Fig. 7.14). These differences can be explained by some uncertainties in the experiments. According to [START_REF] Botia-Vera | Experimental and statistical investigation of canonical problems in sloshing[END_REF], the experiments did not seem to be completely repeatable and symmetric, especially for the shallow oil case where the beam deformations are smaller. Besides, a more precise definition of the material characterization would be needed to better represent the rubber viscoelastic behaviour, which stands for an additional uncertainty. It can also be noted that ∆x F E 0. our 3D solutions is in much better agreement with the experiments that the 2D solution of the literature. 2D simulation leads to a more dynamic flow since the viscous effects on the lateral walls are not considered (see Fig. 7.20).

Here again, good agreements are observed for the free surface deformations (see Fig.

Discussion

A validation of our FSI strategy on 3D cases has been proposed in this chapter. Three different test cases considering free surface flows and deformable bodies have been investigated. For each case, good agreements were obtained between numerical and experimental solutions. Nevertheless, additional experimental data would be needed to validate especially details of the flow evolution. For instance, PIV measurements would be useful to validate the velocity field of our simulations. In addition, the robustness of our coupling method has also been underlined.

Chapter 8

3D aquaplaning problem

The characteristics of our SPH-FE model were detailed in the previous sections. First of all, the FSI coupling strategy has been studied from an energy point of view and computationally improved and validated using 2D test cases. Then the coupling method has been successfully validated using 3D experimental data available in the literature. The purpose of the present chapter is to use this model in order to build a numerical tool able to efficiently model the aquaplaning problem. This challenging problem involves complex 3D fluid-structure interaction. In this chapter, simulations are performed on 3D configurations involving an ideal smooth ground. In this first approach a unique tire design (corresponding to a recent tire) is studied. Note that all following simulations are performed using Abaqus software for the tire.

Model description

The simulation of a tire rolling through a puddle of water initially at rest is proposed (see Fig. 8.1). This simulation aims at analysing how the fluid acts on this complex structure and how the tire deforms due to the presence of water. The eventual persistence of a film of fluid entrapped between the tire and the road is also studied.

In this highly dynamic problem, the use of the SPH method for the fluid description is particularly meaningful. The SPH method naturally considers the three complex interfaces: the free surface, the ground/fluid and the fluid/tire interfaces. It is perfectly suitable for such simulations involving large deformations in both fluid and solid sub-domains. Meshbased method are rather unsuitable for this kind of fluid flow, as the flow physics appears difficult to capture. Indeed, the adaptation between fluid and solid meshes would be very tedious. The SPH method avoids this problem due to the simplification of the fluid-structure interface treatment. Furthermore, mesh-based methods would usually impose a multiphase model with the air phase described, providing additional cells and a dedicated free surface algorithm (VOF or Level-Set). On the contrary, it is possible to model only the water single phase with the SPH method since the free surface conditions are naturally fulfilled with this method. Obviously, it will be needed to investigate up to which extent this single-phase approximation holds in the considered problem.

The puddle is 266 mm long, 180 mm wide and 1 mm high. In the literature higher water heights are usually considered, such as in [START_REF] Kumar | Study of hydroplaning risk on rolling and sliding passenger car[END_REF]. The water height considered here corresponds to the value of the wet braking test imposed by the legislation for tire labelling. It could also be encountered on wet roads during a strong storm. As the wheel enters the puddle, the fluid in front of the tire is entrapped between the tread (see Fig. 8.2) and the ground so that it gets pressurized, leading to the lift generation. The fluid spatial resolution has been set to ∆x SP H = 0.25 mm to correctly capture the pressure field in this area while maintaining reasonable computational costs. The fluid domain is thus discretized with 3.1 millions particles, corresponding to thirty particles per groove width (see Fig. 8.1). It could seem insufficient as it corresponds only to 4 particles in the thickness direction. Nevertheless, this assumption will be validate in the next chapter. The ratio R ∆x SP H has been set to 3.0 to ensure a good compromise between accuracy and needed computational resources. The maximum velocity is expected to be lower than 30 m.s -1 . Therefore, the sound speed is set to 300 m.s -1 to minimize the computational costs while ensuring the weakly-compressible assumption. The fluid simulation parameters are summarized in Tab. 8.1. Note that the complexity of the tire geometry requires the use of the NFM method for the solid boundary conditions in the SPH solver. The tire is modelled using an axi-symmetric tread and carcass (see Fig. 8.2). Among other materials, this complex composite is made of rubber gum, requiring elastic and hyperelastic models. The hyper-elastic model of Mooney-Rivlin is used here. The tread spatial resolution is set to ∆x F E = 2 mm to properly capture the local deformations of the tread induced by the fluid loading, especially in the area where the water is entrapped between the tire and the ground. Such a fine spatial resolution is also needed to correctly handle the solid-solid contact between the tire and the smooth ground. Furthermore, the tire tread Figure 8.3: Free surface shape at t = 10 ms (the tire interface is not represented).

The computational costs are significant due to the complexity of the tire model and the fine fluid spatial resolutions needed in the puddle. Note that the FE solver requires larger CPU times to solve this complex problem than SPH-flow for a single time step. Furthermore, very small time steps (∆t ≈ 10 -6 s) are imposed by the CFL condition of the SPH method, increasing the global computational cost. As a result, the simulation was performed using different time steps in each sub-domain, using the CPS procedure combined with a ratio p = 8. Theoretically, the CPU time in the solid sub-domain (and therefore the total simulation time) is expected to be reduced from a factor 6.4 (Eq. (6.8)) provided that sufficient computing resources are allocated. Here, 120 cores and 16 cores are used respectively for the fluid and the solid sub-domains. In practice, a factor 3.2 has been obtained on the total simulation time. The simulation needed approximately 3 days. Fig. 8.4 shows the time history of the vertical fluid force applied to the tire. For this smooth ground, the support zone (the area where the fluid is mainly taking a part of the solid loading) is located in front of the tread, corresponding to a range of 10% to 20% of the solid loading imposed to the tire. Note that the fluid lift in this configuration thus represents a non-negligible part of the solid loading (337 daN). This magnitude is in agreement with Michelin's expertise for this classical configuration. Furthermore, there is no obvious link between the wet surface and the fluid lift, as the wet surface is monotonically increasing during the simulation (see Fig. 8.5).

From the energetic point of view, an energy transfer occurs from the tire to the puddle of water (see Fig. 8.6). This energy results in fluid particle motions, visible through an increase of the fluid kinetic energy. Water particles are evacuated in front of the tire tread or inside the grooves. Note that the majority of the transmitted energy is dissipated to maintain the stability in the SPH fluid domain, highlighting that the model is far from being converged using this fluid spatial resolution. This was expected since the puddle height is discretized by only 4 particles. Note that following Fig. 8.7 the global system loses energy through the fluid-structure interface. The force conservative formulation described in Chapter 5 is used here. Nevertheless, the loss of energy remains acceptable (less than 2% of E T OT once reached the stationary condition), as discussed in Chapter 5. the tire evacuates the water particles with high blasting velocities. A second portion of water is discharged into the longitudinal grooves. The highest velocities are observed in the jets formed at the beginning of the grooves from the water locally evacuated under the tire. Note that the pressure field is strongly perturbed by the acoustic pressure waves reflected by the tire and the road boundaries which result in global high frequency loading observable in Fig. 8.4. The damping value α s = -0.3 is therefore particularly needed in the HHT scheme for stability purpose (see Section 4.4). Looking carefully at the velocity field, the flow exceeds 30 m.s -1 in some areas (see Fig. 8.9) so that the Mach number is locally higher than 0.1. The validity of our weaklycompressible approach needs therefore to be checked to ensure the relevance of the FSI model. For that purpose, higher values of the sound speed have been tested, C 0 =600 and C 0 =1500 m.s -1 . For this parametric study, the physical durations have also been reduced to 5 and 2 ms respectively to limit the computational costs issues. Fig. 8.10 presents the time history of the resulting vertical fluid forces, showing only small differences for the different sound speed values. These results ensure both the validity of the weaklycompressible approach with C 0 = 300 m.s -1 and the minimization of the number of SPH time steps thanks to this lower sound speed.

According to the fluid pressure field (see Fig. 8.8) the fluid loading is mostly provided by the water in front of the tread. A strong pressure decrease appears in the tire grooves, leading to small fluid forces. The fluid is actually responsible for a loss of the ground reaction force (see Fig. 8.11), which is a key parameter of the grip performances. Note the poor control on this reaction force. Indeed, the model imposes a solid loading of 337 daN and recovers a ground reaction of order 480 daN at the end of the solution, this also without any fluid-structure interaction. The oscillations are also significant, approximatively 200 daN are observed between minimum and maximum values. These force variations are reflected also on the contact surface between the ground and the tire (see Fig. 8.12), which is also a key parameter of the grip performances. This surface is obtained from the Abaqus solver. It corresponds to a summation of the areas where the tire cells intersect the ground cells. Note that both force and surface increase in a non-physical manner.

These differences could be explained by some difficulties in correctly representing the tire model using Abaqus software. Unfortunately, these uncertainties on the solid part prevent the validation of the FSI strategy with experimental data from Michelin's facilities. Indeed, the contact surface between the tread and the ground could be experimentally captured using imaging techniques, as shown in Fig. 8.13. For that purpose, the tire should be correctly modelled. Nevertheless, these uncertainties on the solid modelling do not prevent us from studying the performances of the present FSI coupling and the physic of aquaplaning problems in relative terms. As expected the contact surface is lower with the wet ground. The differences between wet and dry roads can then be studied. It is then possible to analyse (at least partially) the water effects on the global tire grip properties. For example, the total wet surface can be separated into two parts: the wet area in front of the tread and the wet area inside the grooves (see Fig. 8.9). After t = 5 ms, the wet surface area in front of the tire remains globally constant, as shown in Fig. 8.14. A rather steady state is quickly obtained with regard to the evolution of this area. Conversely, no steady state is reached in the grooves yet. Water jets need longer times to propagate inside the tire structures. The simulation is extended to 10 ms in order to let the flow to fully develop, but this is still not sufficient. Nevertheless, it does not influence the fluid lift since the contribution of the grooves is limited here to the fluid evacuation. For wet and dry grounds, averaged values of the vertical solid force and of the solid-solid contact surface are estimated for t ≥ 5 ms. For t ≤ 5 ms, the flow is not established yet (looking at the wet area in front of the tire). It is then possible to estimate the variations between wet and dry grounds (see Tab. 8.2). These results are in global agreement with the experimental behaviour and Michelin's expertise for such aquaplaning configurations. It is important to emphasize that some data is almost impossible to obtain experimentally with the technologies currently available. This FSI coupling is therefore very useful to understand the aquaplaning. The SPH-FE coupling can thus complete local or global experimental quantities for tire design purpose. For example, the less the water is evacuated into the grooves, the more it moves towards the water bead in front of the contact area. The fluid is then intensely pressurised. As a result, the tire is locally deformed by these strong fluid stress. We propose here to study the local deformations at three rib cross-sections (of 4 mm width), as described in Fig. 8.15. Fig. 8.16 presents the results at t = 2.5 ms (x and z correspond to the global coordinates presented in Fig. 8.1). The puddle of water tends to locally press the front of the tire. On a dry ground, the tire profile is expected to be the same at the front and at the back of the contact area. The deformations are larger in the lateral cross-sections than in the central one, which contradicts the experimental feedback. Nevertheless, this difference can be easily explained by the difficulties to correctly represent the tire problem with Abaqus, especially the mechanical behaviour of each tread component.

Averaged

Note also that the loss of contact surface is observed exclusively at the front (the back remains unchanged). Fig. 8.17 displays some snapshots of the contact pressure field due to the ground on the tire. Pressures have been limited to 3.10 5 Pa to make the observation easier and the tire moves from the right to the left. Looking carefully at these snapshots, the pressure distribution is only modified at the front (left side). Indeed, in this area the water tends to sustent the tire from the ground, reducing the contact pressure. This local elevation of the tire located in front of the contact area is responsible for a modification of the vertical displacement of the wheel center, as shown is Fig. 8.18. Note that the dry displacement is already decreasing, due to the difficulties to properly modelled the tire structure in Abaqus. Afterwards, we realized that the tire initialization should contain a acceleration phase on a dry ground to remove any potential transient effects. However, this step seems impossible to perform if the tire tread is only represented on 70 degrees. Nevertheless, a whole wheel simulation with fine tire spatial resolution is unusable for computational time reasons. At t = 10 ms, the wheel center has moved of about 0.4 mm due to the fluid forces. This relative displacement remains small since a low water height is considered. This data is interesting for the validation of such aquaplaning simulations. Indeed, this displacement can be measured experimentally at Michelin's facilities. These simulations provide therefore an order of magnitude of the displacement, which helps purposing similar experiments.

Adaptive particle refinement

The number of particles involved in the previous simulations is too large (at least few millions) to ensure lower CPU times in SPH than in FE (with reasonable numbers of CPU cores):

p t SP H > t F E . (8.1)

However, this condition is necessary to approach the maximum theoretical gain. The fluid spatial resolution ∆x SP H = 0.25 mm is needed in the near-body area to correctly describe the pressure field. Nevertheless, the particles located upstream the puddle do not need such a spatial resolution. In the aquaplaning problem, the area of interest is not fixed since the tire is moving into the puddle. As a result, the particle refinement method can be used to limit the number of particles involved according to the configuration presented in Fig. 8.19 (∆x SP H = 0.50 mm in the non-refined area) while preserving the desired accuracy. The refinement box should therefore follow the pressurized area in front of the tire. The aim is to maintain the SPH computational time of p time steps lower than the FE one for a single time step (p t SP H < t F E ).

The resulting number of particles is divided by approximatively three. Note that this method is increasingly beneficial with the puddle length. The same tire model is adopted here and the fluid parameters are summarized in Tab. 8.3. Note that the APR technique needs the particle shifting (described in Section 2.7), explaining the differences with the previous results. Consequently, a fully refined simulation is also performed to serve as reference case. Nevertheless the fluid flow is strongly modified outside the refinement box due to a loss of accuracy in the fluid description. Fig. 8.25 displays some snapshots of the computed velocity field with and without APR, showing similar fields inside the refined area. Note that the free surface outside the APR box is altered by the coarse fluid spatial resolution. Consequently, it is not possible to obtain fully satisfactory results concerning the wet surface area or the flow rate inside the grooves. For that purpose, several refinement boxes are needed. Otherwise, the refinement box would have to be extended to the corresponding areas. However it does not have any influence on the global fluid-structure interaction (prediction of fluid loads and corresponding deformations) since it corresponds to low pressurized area. Concerning the global computational time, the practical gain has been up to a factor 3.7 compared to a simulation using p = 1 and without APR box. It runs during 2.6 days using respectively 120 and 16 cores for the fluid and the solid solvers. Otherwise the practical gains were limited to 3.2 (i.e. 2.9 days) in the fully refined configurations. The time saving of the APR method is not very significant in this configuration, a decrease of 15%. Note that the performances of this method is directly linked to the number of fluid particles per processor. A low number has been used here using 120 cores with SPH-flow. Nevertheless, the APR method has the additional advantage of reducing the volume of fluid data. The total volume is approximately divided by 3. This is not negligible for data post-processing purpose with several millions of SPH particles, especially in an industrial context.

Discussion

The development of a complete tire model has been performed. This model is quite representative: pre-inflation, solid loading applied to the center of the wheel, initial contact setting, rotation and translation speeds imposed to the rim flank. The capability of the current FSI strategy to deal with this multi-dimensional problem involving complex interfaces has been demonstrated. The present FSI strategy naturally considers the complex interfaces without any particular issue. The first simulations have provided coherent results for this smooth ground configuration. It has also been underlined that several parameters should be studied (fluid lift, wet surfaces, flow rates, tire deformations, etc.), helping in understanding the mechanisms at the origin of aquaplaning. The aim of this thesis was to develop an aquaplaning model to improve the understanding of its generating mechanisms. A weak partitioned coupling based on the SPH and FE methods has been chosen to model this challenging fluid-structure interaction problem. This strategy is built on appropriate exchanges between two distinct software, allowing to use the best advantages of both methods for each sub-domain.

Firstly, an energetic study of the coupling strategy has been performed using 2D test cases. A particular attention has been paid to the interface coupling energy. This allowed us to validate the FSI coupling from an energetical point of view. It also permitted the improvement of its robustness and accuracy properties using new fluid loading on the deformable structure.

Various tools were proposed and tested to reduce the computational costs of such simulations. Industrial applications require reasonable computational times to be efficient. The interest to use different time steps in each sub-domain has been demonstrated, mainly through the use of two coupling schemes. Furthermore, the particle refinement showed some good capabilities to reduce the coupling computational time without altering its accuracy and stability properties.

The relevance of our 3D SPH-FE model has been investigated. The coupling strategy has been therefore benchmarked using various experimental data. Several test cases involving different flow dynamics and material deformations have been successfully simulated. Good agreements with available experimental data have been found in each configuration.

The feasibility of 3D aquaplaning simulations on a ideal smooth ground using the SPH-FE coupling strategy has been analysed. In each case, the numerical results show good agreement with experimental behaviours observed at Michelin's facilities. This strategy provides a large variety of outputs helping in the understanding of the aquaplaning generating mechanics.

CHAPTER 11. CONCLUSIONS AND PERSPECTIVES

Perspectives

The work covered in this thesis was an introduction to a large range of possible researches on the aquaplaning phenomenon. It proves the feasibility of such simulations using a SPH-FE coupling strategy. Nevertheless, there are many prospects for improving the FSI model since we have been confronted to various challenges during this work.

The major perspective would be to improve the model on the solid part. Indeed Abaqus software has demonstrated some difficulties to correctly represent the tire mechanics. Another coupling is therefore needed between SPH-flow and Michelin in-house software dedicated to the FE simulation. Once done, the proposed model will have to be applied to a conventional (reference) tire to assess the model accuracy using the experimental database from Michelin's facilities. Some points deserve further interrogations: description of friction, larger rotations of the tire, more complicated groove designs such as transversal of non-symmetric sculptures such as in OKANO andKOISHI (2001) or NAKAJIMA et al. (2000), materials, etc.

Moreover for different reasons, some questions have been left unanswered but could be pursued:

• The possible influence of the ventilation has to be considered. Indeed, some air pockets may persist inside the tire grooves, preventing the water from being correctly evacuated. The grip properties could thus be modified by considering the air.

• Another perspective is to consider the possible influence of turbulence and/or cavitation effects into the fluid flow, even though the SPH method does not seem particularly suitable for such modelling.

Few attempts have been made to the study the ventilation and cavitation effects but without convincing results. That is why they have not been presented in this manuscript.

In addition, note that all specific improvements of the SPH or FE methods could be beneficial to the coupling strategy. For an example, improvements providing more regular pressure fields would be particularly useful for stability issues.

Finally, an extension of the present work is obviously to pursue the validation of the SPH-FE coupling strategy on more complex fluid-structure interaction problems, involving non-linearity, large displacements, etc.

F +S E 0 F

 0 Time history of ratio between E Interf ace and E T OT for different stabilization coefficients. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.12 Time history of ratio between E Interf ace and E V for different stabilization coefficients, ∆x SP H = 1mm. . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.13 Time history of the vertical force on the beam for different stabilization coefficients, ∆x SP H = 1mm. . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.14 Pressure field comparisons between different stabilization coefficients (rows) for different instants (columns), ∆x SP H = 1 mm. . . . . . . . . . . . . . . . 5.15 Dam-break configuration. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.16 Strain-stress curve for the rubber used. . . . . . . . . . . . . . . . . . . . . . 5.17 Time history of the horizontal (top) and vertical (bottom) displacements of the gate tip for different fluid resolutions compared to experiments, Newmark scheme with α s = -0.3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.18 Time history of solid total energy for different fluid spatial resolutions, Newmark scheme with α s = -0.3. . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.19 Time history of E M (dashed lines), E K (dotted lines) and E P (lines) for different fluid spatial resolutions, Newmark scheme with α s = -0.3. . . . . . . . 5.20 Time history of ratio between E f luid/body (dashed lines), E body/f luid (lines) for different fluid spatial resolutions, Newmark coefficient α s = -0.3. . . . . . . . 5.21 Time history of E F +S -E 0 +S for different fluid spatial resolutions, Newmark scheme with α s = -0.3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.22 Time history of ratio between E Interf ace and E T OT for different stabilization coefficients, ∆x SP H = 0.25mm. . . . . . . . . . . . . . . . . . . . . . . . . . 5.23 Time history of the horizontal (top curve) and vertical (bottom curve) displacements of the gate tip for different stabilization coefficients, ∆x SP H = 0.25mm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.24 Pressure field comparisons between different stabilization coefficients (rows) at different instants (columns), ∆x SP H = 0.25 mm. . . . . . . . . . . . . . . 5.25 Time history of ratio between E Interf ace and E V for different stabilization coefficients, ∆x SP H = 0.25 mm. . . . . . . . . . . . . . . . . . . . . . . . . . 5.26 Time history of E V for different fluid spatial resolutions. . . . . . . . . . . . . 5.27 Time history of ratios between E Interf ace and E T OT (left) and between E Interf ace and E V (right) for different stabilization coefficients, force conservative formulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.28 Time history of ratios between E Interf ace and E T OT (left) and between E Interf ace and E V (right) for different stabilization coefficients, energy conservative formulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . LIST OF FIGURES 5.29 Time history of vertical force on the deformable beam for different load formulations and boundary conditions, HHT scheme with α s = -0.1. . . . . 98 5.30 Pressure field comparison for different load formulations (columns) at different instants (lines), HHT scheme with α s =-0.05. . . . . . . . . . . . . . . . . 100 5.31 Time history of ratios between E Interf ace and E V for different stabilization coefficients. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 5.32 Time history of ratios between E Interf ace and E T OT for different stabilization coefficients. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 5.33 Time history of horizontal (top) and vertical (bottom) displacements of the gate tip for the different formulations proposed, HHT scheme with α s = -0.1. 102 5.34 Time history of horizontal (top) and vertical (bottom) displacements of the gate tip for the different formulations proposed, Newmark scheme with α s = -0.3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 5.35 Time history of horizontal force on the elastic gate for the different formulations proposed, Newmark scheme with α s = -0.3. . . . . . . . . . . . . . . . 102 5.36 Pressure field comparison for different load formulations (columns) at different instants (lines), HHT scheme with α s =-0.05. . . . . . . . . . . . . . . . . 103 6.1 Coupling algorithms with different time steps. . . . . . . . . . . . . . . . . . 107 6.2 Time history of vertical force on the deformable beam for different p ratios and coupling algorithms. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 6.3 Time history of the midpoint relative displacement on the deformable beam for different p ratios and coupling algorithms. . . . . . . . . . . . . . . . . . 109 6.4 Pressure field comparisons for different p ratios at different instants, CPS procedure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 6.5 Pressure field comparisons for different p ratios at different instants, CPS procedure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 6.6 Pressure field comparison at t = 3 ms for different p ratios, CPS procedure. . 112 6.7 Pressure field comparison at t = 3 ms with and without linear extrapolation Eq. (6.1) for p = 16, CPS procedure. . . . . . . . . . . . . . . . . . . . . . . 113 6.8 Time histories of the vertical force (left) and the midpoint relative displacement (right) on the deformable beam for different p ratios without extrapolation Eq. (6.1), CPS procedure. . . . . . . . . . . . . . . . . . . . . . . . . . 113 6.9 Time history of the horizontal force on the elastic gate for different p ratios and coupling algorithms. . . . . . . . . . . . . . . . . . . . . . . . .

7. 1

 1 3D dam-break configuration. . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2 Time history of the horizontal (top) and vertical (bottom) displacements of the gate tip for different FSI models compared to experiments and to LI et al. (2015). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.3 Water height evolutions at locations x = 0.4 m (left) and x = 0.45 m (right) for different FSI models compared to experiments. . . . . . . . . . . . . . . . 7.4 Fluid pressure field at different instants (rows) for different FSI models compared to corresponding snapshots from Antoci et al. experiments. . . . . . . 7.5 Dam-break configuration with a single APR box (red). . . . . . . . . . . . . 7.6 Time history of the horizontal (top) and vertical (bottom) displacements of the gate tip for different FE solvers with APR compared to experiments. . . 7.7 Water height evolutions at locations x = 0.4 m (left) and x = 0.45 m (right) for different FE solvers with APR compared to experiments. . . . . . . . . . 7.8 Fluid pressure field at different instants (rows) for different FE solvers with APR compared to corresponding snapshots from Antoci et al. experiments, Y = 0.05 m. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.9 Two initial configurations of the FSI sloshing problem. . . . . . . . . . . . . 7.10 Rolling tank configuration, BOTIA-VERA (2015). . . . . . . . . . . . . . . . 7.11 Gaps between the clamped beam and the tank. . . . . . . . . . . . . . . . . 7.12 Clamped beam in mid depth oil configuration. . . . . . . . . . . . . . . . . . 7.13 Time history of the forced roll motion. . . . . . . . . . . . . . . . . . . . . . 7.14 Time history of the local x-displacement of the the beam extremity. . . . . . 7.15 Current simulation results (left) compared to the experiments (right) at different instants (rows). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138 7.16 Pressure field comparisons for 2D (left) and 3D (right) simulations at different instants. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139 7.17 Clamped beam in shallow oil configuration. . . . . . . . . . . . . . . . . . . . 140 7.18 Time history of the tank inclination. . . . . . . . . . . . . . . . . . . . . . . 141 7.19 Time history of the local X-displacement of the the beam extremity. . . . . . 141 7.20 Current simulation results (left) compared to the experiments (right) at t = 1.68 s for 2D and 3D simulations. . . . . . . . . . . . . . . . . . . . . . . . . 141 7.21 Current simulation results (left) compared to the experiments (right) at different instants (rows). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142 8.1 Aquaplaning configuration with the smooth ground, the puddle initially at rest and the tire interface. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146 8.2 Tire tread (red) and carcass (green) used in the simulation. . . . . . . . . . . 147 8.3 Free surface shape at t = 10 ms (the tire interface is not represented). . . . . 148 8.4 Time history of the vertical fluid force acting on the tire, smooth road. . . . 149 8.5 Time history of the tire total wet surface, smooth road. . . . . . . . . . . . . 149 8.6 Time history of E V , E K and E T OT in the fluid, smooth road. . . . . . . . . . . 149 8.7 Time history of the energy ratio between E Interf ace and E T OT , smooth road. . 149 8.8 Velocity (left) and pressure (right) fields at different instants, smooth road. . 150 8.9 Zoom of the velocity field at t = 2 ms, smooth road. . . . . . . . . . . . . . . 151 8.10 Time history of the vertical fluid force acting on the tire for different sound speed, smooth road. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152 8.11 Time history of the vertical solid force on tire for wet and dry smooth roads. 152 8.12 Time history of the solid-solid contact surface for wet and dry smooth roads. 152 8.13 Aquaplaning experimental tests, smooth road. . . . . . . . . . . . . . . . . . 153 8.14 Time history of the tire wet surfaces, smooth road. . . . . . . . . . . . . . . 153 8.15 The rib cross-sections. . . . . . . . . . . . . . . . . . . . . . . . . .

A∇

  The dissipative coefficient of the δ-SPH scheme α AV The artificial viscosity coefficient α s Stability parameter of the solid time integration scheme β The kernel normalization constant β l The minmond limiter coefficient β s Stability parameter of the solid time integration scheme ∆t The time step ∆x The distance between two particles ∆x c The child particle spacing for APR ∆x p The parent particle spacing for APR δ The Dirac function δ δ The diffusive coefficient of the δ-SPH scheme The The matrix of the stress-strain constitutive relation ∂D The contour of the kernel support ∂Ω The boundary of the domain Ω ∂Ω s The solid domain boundary ∂Ω t The solid domain for the strain constraints ∂Ω v The solid domain for the velocity constraints RL W The renormalization kernel gradient φ The vector of conservative variables F υ The viscous flux tensor F c The convective flux tensor F body/f luid The force applied by the deformable body on the fluid F f luid/body The force applied by the fluid on the deformable body F ext The generalized force vector Monotone Upstream-Centered Scheme for Conservative Laws N The number of particles in the support kernel N F E The number of Finite Element N g The number of ghost particles in the kernel support N p The number of wall panels inside the kernel support The tire rigidity relative to the water height S p The speedup SP H Smoothed Particle Hydrodynamics t F E The CPU time of the FE solver per time step t SP H The CPU time of the SPH solver per time step

Figure 1 . 1 :

 11 Figure 1.1: Aquaplaning transition zones.

  .14) where D stands for the compact support of the kernel function W. The size of D depends on h, the smoothing length. For the present manuscript, the radius R of the compact support W (proportional to h) is introduced. A dependency is thus created from the considered point and its immediate neighbourhood, Fig. 2.1.

Figure 2 . 1 :

 21 Figure 2.1: Schematic representation of the kernel and its compact support.

CHAPTER 2 .

 2 Figure 2.2: Arbitrary particle distribution in a kernel support.

Figure 2 .

 2 Figure 2.3: Fluid domain and its various boundary conditions.

Figure 2 . 4 :

 24 Figure 2.4: Ghost particle method.

Figure 2 . 5 :

 25 Figure 2.5: Normal flux method.

  2.6).

Figure 2 .

 2 Figure 2.6: Free surface configurations.

Figure 2 .

 2 Figure 2.7: Particle refinement process : parent particles (in red) are split into four child particles (in blue).

CHAPTER 2 .

 2 Figure 2.8: The particle refinement process.

  T b dS.The virtual displacement is chosen as null on the boundary ∂Ω v :

Figure 4 . 1 :

 41 Figure 4.1: NFM conditions for SPH-FE coupled simulations.

Figure 4 . 2 :

 42 Figure 4.2: Averaged pressure calculation on a body panel k.

Figure 4 .

 4 Figure 4.3: Coupling algorithms.

Figure 4 . 4 :

 44 Figure 4.4: Hydrostatic water column on an elastic plate.

L

Figure 4 . 5 :

 45 Figure 4.5: Vertical displacement time history of the beam midpoint for different CFL values, Newmark scheme with α s = 0.

Figure 4 .

 4 Figure 4.6: Vertical displacement time history of the beam midpoint for different stabilizing coefficients, HHT scheme.

Figure 4

 4 Figure 4.7: Pressure field comparison for different values of α s in the HHT scheme (rows) at different instants (columns) using the CPS algorithm, only the perturbation on pressure is represented (i.e. the hydrostatic component is subtracted).

  [START_REF] Antuono | Energy balance in the δ-SPH scheme[END_REF] proposed a study related to energy conservation in the SPH method, in the absence of solid bodies and for the δ-SPH scheme. Their work analyses CHAPTER 5. ENERGY CONSIDERATIONS ON THE SPH-FE COUPLING STRATEGY the contribution of this term from the continuity equation, highlighting that the diffusive term is generally small and linked to the excitation of spurious high-frequency acoustic waves. Then CERCOS-PITA et al. (2017) performed this investigation in presence of rigid body interactions, highlighting that imposing boundary conditions using ghost particles introduces some extra-terms in the energy equations. These extra-terms are actually related to the fluid-structure interactions. Two different boundary techniques are investigated here, the ghost particles method (Sec. 2.5.1.1) and the Normal Flux Method (Sec.2.5.1.2).

Figure 5 . 1 :

 51 Figure 5.1: Energy balance over the global domain.

5. 1 .

 1 ENERGY CONSIDERATIONS OVER THE GLOBAL DOMAIN dE f luid/body dt = k∈B P k S k n k . v k , (5.1) where B is the set of wet faces of surface S k , normal n k , local velocity v k and P k is the fluid pressure applied on it (Fig. 4.2).

Figure 5 . 2 :

 52 Figure 5.2: Fluid-structure interface from the solid (left) and fluid (right) point of views.

Figure 5 .

 5 Figure 5.3: Deformable beam impact configuration.

  5.5).

Figure 5 . 4 :

 54 Figure 5.4: Time history of the midpoint relative displacement for different spatial resolutions in the fluid, Newmark scheme with α s = -0.3.

Figure 5 . 5 :

 55 Figure 5.5: Time history of vertical force on the deformable beam for different spatial resolutions in the fluid, Newmark scheme with α s = -0.3.

Figure 5 .

 5 Figure 5.6: Time history of the pressure probes on the deformable beam for different spatial resolutions in the fluid, Newmark scheme with α s = -0.3.

Figure 5 . 7 :

 57 Figure 5.7: Time history of E M (lines) and E body/f luid (dashed lines) of the water for different spatial resolutions in the fluid, Newmark scheme with α s = -0.3.

Figure 5 . 8 :

 58 Figure 5.8: Time history of E V (dashed lines) and E C (lines) of the water for different spatial resolutions in the fluid, Newmark scheme with α s = -0.3.

Figure 5 .

 5 Figure 5.9: Time history of E f luid/body (lines) and E body/f luid (dashed lines) for different spatial resolutions in the fluid, Newmark scheme with α s = -0.3.

Figure 5 .

 5 Figure 5.10: Time history of ∆E F +S (lines) and E Outside (dashed lines) for different spatial resolutions in the fluid, Newmark scheme with α s = -0.3.

Figure 5 .

 5 Figure 5.11: Time history of ratio between E Interf ace and E T OT for different stabilization coefficients.

Figure 5 . 12 :

 512 Figure 5.12: Time history of ratio between E Interf ace and E V for different stabilization coefficients, ∆x SP H = 1mm.

Fig. 5 .

 5 Fig. 5.13 compares the vertical force time histories for HHT and Newmark schemes for different values of the dissipation parameter α s . High frequency oscillations are observed for low dissipation values in the HHT scheme, although the energy evolution matches with the coupling stability criterion described in FOUREY et al. (2017).Fig.5.14 confirms the above statements, by presenting snapshots of the computed solution at various instants with different FE time integrator schemes. These snapshots are plotted with pressure and Von Mises stress contours respectively for the fluid and the structure. Strong non physical perturbations are observed in the pressure field for HHT with α s = -0.05 (especially after t=2ms), leading to strong vertical oscillations as visible in Fig.5.13. In this configuration, increasing the stabilization parameter α s helps in maintaining a stable coupling and reducing the energy errors (Fig.5.12).

5. 4 .

 4 Figure 5.13: Time history of the vertical force on the beam for different stabilization coefficients, ∆x SP H = 1mm.

CHAPTER 5 .Figure 5 .

 55 Figure 5.14: Pressure field comparisons between different stabilization coefficients (rows) for different instants (columns), ∆x SP H = 1 mm.

Fig. 5 .

 5 Fig.5.22 shows the ratio between E Interf ace and E T OT with different solid time integration

Figure 5 .

 5 Figure 5.15: Dam-break configuration.

Figure 5 .

 5 Figure 5.16: Strain-stress curve for the rubber used.

Figure 5 .

 5 Figure 5.17: Time history of the horizontal (top) and vertical (bottom) displacements of the gate tip for different fluid resolutions compared to experiments, Newmark scheme with α s = -0.3.

Figure 5 .

 5 Figure 5.18: Time history of solid total energy for different fluid spatial resolutions, Newmark scheme with α s = -0.3.

Figure 5 .

 5 Figure 5.19: Time history of E M (dashed lines), E K (dotted lines) and E P (lines) for different fluid spatial resolutions, Newmark scheme with α s = -0.3.

Figure 5 .

 5 Figure 5.20: Time history of ratio between E f luid/body (dashed lines), E body/f luid (lines) for different fluid spatial resolutions, Newmark coefficient α s = -0.3.

Figure 5

 5 Figure 5.21: Time history ofE F +S -E 0 F +S

Figure 5 .

 5 Figure 5.23: Time history of the horizontal (top curve) and vertical (bottom curve) displacements of the gate tip for different stabilization coefficients, ∆x SP H = 0.25mm.

Figure 5 .

 5 Figure 5.24: Pressure field comparisons between different stabilization coefficients (rows) at different instants (columns), ∆x SP H = 0.25 mm.

Figure 5 .

 5 Figure 5.25: Time history of ratio between E Interf ace and E V for different stabilization coefficients, ∆x SP H = 0.25 mm.

Figure 5 .

 5 Figure 5.26: Time history of E V for different fluid spatial resolutions.

Figure 5 .Figure 5 .

 55 Figure 5.30: Pressure field comparison for different load formulations (columns) at different instants (lines), HHT scheme with α s =-0.05.

Figure 5 .

 5 Figure 5.32: Time history of ratios between E Interf ace and E T OT for different stabilization coefficients.

Figure 5 .

 5 Figure 5.33: Time history of horizontal (top) and vertical (bottom) displacements of the gate tip for the different formulations proposed, HHT scheme with α s = -0.1.

Figure 5 .

 5 Figure 5.34: Time history of horizontal (top) and vertical (bottom) displacements of the gate tip for the different formulations proposed, Newmark scheme with α s = -0.3.

Figure 5 .Figure 5 .

 55 Figure 5.35: Time history of horizontal force on the elastic gate for the different formulations proposed, Newmark scheme with α s = -0.3.

6. 1

 1 SPH-FE coupling using different time steps in each sub-domain 6.1.1 Coupling algorithm details The optimization proposed here consists in authorizing several fluid time steps within each solid one. Recently NUNEZ RAMIREZ et al. (2016) proposed such a coupling strategy able to integrate each sub-domain with different time steps by adapting the energy-conserving procedure proposed by LI et al. (2015). However, partitioned approaches (CPS and CSS) are still considered in this work. As a consequence, the Conventional Parallel Staggered and Conventional Sequential Staggered procedures should be adapted. Here we propose to introduce a ratio p between the numbers of fluid and solid time steps (see Fig. 6.1).

6. 1 .Figure 6 . 2 :Figure 6 . 5 :

 16265 Figure 6.2: Time history of vertical force on the deformable beam for different p ratios and coupling algorithms.

Figure 6

 6 Figure 6.6: Pressure field comparison at t = 3 ms for different p ratios, CPS procedure.

  Figure 6.9: Time history of the horizontal force on the elastic gate for different p ratios and coupling algorithms.

  Figure 6.13: Pressure field comparisons for different p ratios at different instants, CSS procedure.

Figure 6 .

 6 Figure 6.14: Deformable beam impact configuration with APR boxes (red).

Figure 6 .

 6 Figure 6.15: Time history of vertical force on the deformable beam for different numbers of refinement fixed boxes.

Figure 6 .

 6 Figure 6.16: Time history of the midpoint relative displacement on the deformable beam for different numbers of refinement fixed boxes.

Figure 6 .

 6 Figure 6.17: Time history of vertical force on the deformable beam for different numbers of refinement moving boxes.

Figure 6 .

 6 Figure 6.18: Time history of the midpoint relative displacement on the deformable beam for different numbers of refinement moving boxes.

Fig. 6 .Figure 6 Figure 6

 666 Fig.6.19 shows some snapshots of the computed solution for different APR configurations, highlighting the efficiency of the APR method to reduce the number of particles involved in a simulation without altering the accuracy and stability properties. Further-

Figure 6 .

 6 Figure 6.21: Dam-break configuration with APR boxes (red).

Figure 6 .

 6 Figure 6.22: Time history of horizontal (top) and vertical (bottom) displacements of the gate tip for different numbers of refinement boxes.

Figure 6 .

 6 Figure 6.23: Time history of the horizontal force on the elastic gate for different numbers of refinement boxes.

Figure 6 .

 6 Figure 6.24: Time history of horizontal force on the elastic gate for different numbers of refinement boxes and p ratio.

Figure 6 .

 6 Figure 6.25: Time history of horizontal force on the elastic gate for different p ratios and without refinement boxes.

  Figure 7.1: 3D dam-break configuration.

Figure 7 . 2 :

 72 Figure 7.2: Time history of the horizontal (top) and vertical (bottom) displacements of the gate tip for different FSI models compared to experiments and to LI et al. (2015).

Figure 7 . 3 :

 73 Figure 7.3: Water height evolutions at locations x = 0.4 m (left) and x = 0.45 m (right) for different FSI models compared to experiments.

Figure 7 . 4 :

 74 Figure 7.4: Fluid pressure field at different instants (rows) for different FSI models compared to corresponding snapshots from Antoci et al. experiments.

Figure 7 . 5 :

 75 Figure 7.5: Dam-break configuration with a single APR box (red).

Fig. 7

 7 Fig. 7.6 shows the time history of the gate tip displacement for both Abaqus and Code Aster FE solvers. The FSI behaviour is correctly modelled for each configuration, whatever the FE solver or the refinement level. Note that the displacement is just slightly modified. On the contrary, the water height evolutions and the fluid pressure field are very similar in all cases (see Fig.7.7 and Fig.7.8). These results ensure both validity of the APR method on 3D FSI problems and the coupling robustness regarding the different FE solvers used.

Figure 7 . 6 :

 76 Figure 7.6: Time history of the horizontal (top) and vertical (bottom) displacements of the gate tip for different FE solvers with APR compared to experiments.

Figure 7 . 7 :Figure 7

 777 Figure 7.7: Water height evolutions at locations x = 0.4 m (left) and x = 0.45 m (right) for different FE solvers with APR compared to experiments.

Figure 7 . 9 :

 79 Figure 7.9: Two initial configurations of the FSI sloshing problem.

7. 2 .

 2 Figure 7.10: Rolling tank configuration, BOTIA-VERA (2015).

Figure 7 .

 7 Figure 7.11: Gaps between the clamped beam and the tank.

Figure 7 . 12 :

 712 Figure 7.12: Clamped beam in mid depth oil configuration.

7. 2 .

 2 SLOSHING TANK INTERACTIONS WITH AN ELASTIC BEAMIn this case, the amplitude of the roll motion is about 4 degrees (Fig.7.13). The clamped beam therefore presents large deformations. Fig.7.14 shows the time history of the local x-displacement of the beam extremity. Numerical results are in good agreements with the experimental solution from BOTIA-VERA (2015) and the numerical results from the literature, both in phase and in amplitude.

Figure 7 .

 7 Figure 7.13: Time history of the forced roll motion.

Figure 7 .

 7 Figure 7.14: Time history of the local xdisplacement of the the beam extremity.

Fig. 7 .

 7 Fig.7.15 displays the fluid pressure field and compares the current results to the experiments at different instants. A good agreement is observed for the free surface deformations. This problem seems correctly predicted by the proposed SPH-FE coupling. However, finer experimental data would be needed detailed validation of the numerical model, such as PIV measurements for the velocity field and pressure sensors. The use of a 3D model is crucial here. The fluid flow is much more dynamic in 2D simulations since the viscous effects on the lateral walls are not considered (see Fig.7.16).

  Figure 7.17: Clamped beam in shallow oil configuration.

Figure 7 .

 7 Figure 7.18: Time history of the tank inclination.

Figure 7 .

 7 Figure 7.19: Time history of the local Xdisplacement of the the beam extremity.

Figure 8 . 1 :

 81 Figure 8.1: Aquaplaning configuration with the smooth ground, the puddle initially at rest and the tire interface.

  Figure 8.4: Time history of the vertical fluid force acting on the tire, smooth road.

Figure 8 . 5 :

 85 Figure 8.5: Time history of the tire total wet surface, smooth road.

Figure 8 . 6 :

 86 Figure 8.6: Time history of E V , E K and E T OT in the fluid, smooth road.

Figure 8 . 7 :

 87 Figure 8.7: Time history of the energy ratio between E Interf ace and E T OT , smooth road.

  Figure 8.8: Velocity (left) and pressure (right) fields at different instants, smooth road.

Figure 8

 8 Figure 8.9: Zoom of the velocity field at t = 2 ms, smooth road.

Figure 8 .

 8 Figure 8.10: Time history of the vertical fluid force acting on the tire for different sound speed, smooth road.

Figure 8 .

 8 Figure 8.11: Time history of the vertical solid force on tire for wet and dry smooth roads.

Figure 8 . 12 :

 812 Figure 8.12: Time history of the solid-solid contact surface for wet and dry smooth roads.

Figure 8 .

 8 Figure 8.13: Aquaplaning experimental tests, smooth road.

Figure 8 .

 8 Figure 8.14: Time history of the tire wet surfaces, smooth road.

Figure 8 .

 8 Figure 8.15: The rib cross-sections.

Figure 8 .

 8 Figure 8.18: Time history of the wheel center vertical displacement for wet and dry smooth roads.

Figure 8 .

 8 Figure 8.19: Aquaplaning configuration with a single APR box.

Figure 8 .

 8 Figure 8.20: Time history of the vertical fluid force with and without APR, smooth road.

Figure 8 .

 8 Figure 8.21: Time history of wheel center vertical displacement with and without APR, smooth road.

Figure 8 .

 8 Figure 8.22: Time history of the vertical solid force on tire with and without APR, smooth road.

Figure 8 .

 8 Figure 8.23: Time history of the solidsolid contact surface with and without APR, smooth road.

Figure 8

 8 Figure 8.24: Tire profiles at the left (top), central (middle) and right (bottom) rib crosssections at t = 2.5 ms with and without APR, wet smooth road.

Figure 8

 8 Figure 8.25: Velocity field at different instants with (left) and without (right) APR, smooth road.

A. 2 .

 2 Figure A.4: Time history of the ratio between E Interf ace and E T OT for different stabilization coefficients, energy conservative formulation.

Figure A. 5 :

 5 Figure A.5: Time history of the ratio between E Interf ace and E V for different stabilization coefficients, energy conservative formulation.

Figure A. 6 :

 6 Figure A.6: Time history of the ratio between E Interf ace and E V for different boundary conditions and fluid loading, HHT scheme with α s = -0.1.

Figure A. 7 :Figure B. 2 :

 72 Figure A.7: Time history of the atio between E Interf ace and E T OT for different boundary conditions and fluid loading, HHT scheme with α s = -0.1.

  . . . . . 114 6.10 Pressure field comparisons for different p ratios at different instants, CPS procedure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 6.11 Time history of horizontal (top) and vertical (bottom) displacements of the gate tip for different p ratios and coupling algorithms. . . . . . . . . . . . . . 116 6.12 Pressure field comparisons for different p ratios at different instants, CSS Time history of vertical force on the deformable beam for different numbers of refinement fixed boxes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.16 Time history of the midpoint relative displacement on the deformable beam for different numbers of refinement fixed boxes. . . . . . . . . . . . . . . . . 6.17 Time history of vertical force on the deformable beam for different numbers of refinement moving boxes. . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.18 Time history of the midpoint relative displacement on the deformable beam for different numbers of refinement moving boxes. . . . . . . . . . . . . . . . 6.19 Pressure field comparisons for different APR configurations. . . . . . . . . . 6.20 Pressure field comparisons for different APR configurations. . . . . . . . . . 6.21 Dam-break configuration with APR boxes (red). . . . . . . . . . . . . . . . . 6.22 Time history of horizontal (top) and vertical (bottom) displacements of the gate tip for different numbers of refinement boxes. . . . . . . . . . . . . . . . 6.23 Time history of the horizontal force on the elastic gate for different numbers of refinement boxes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.24 Time history of horizontal force on the elastic gate for different numbers of refinement boxes and p ratio. . . . . . . . . . . . . . . . . . . . . . . . . . . 6.25 Time history of horizontal force on the elastic gate for different p ratios and without refinement boxes. .

	LIST OF FIGURES
	6.15

procedure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 6.13 Pressure field comparisons for different p ratios at different instants, CSS procedure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 6.14 Deformable beam impact configuration with APR boxes (red). . . . . . . . . 119 . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.26 Pressure field comparisons at different instants (columns) for different number of refinement fixed boxes and different ratios p. . . . . . . . . . . . . . .

  . . . . . 155 8.16 Tire profiles at the left (top), central (middle) and right (bottom) rib crosssections at t = 2.5 ms, wet smooth road. . . . . . . . . . . . . . . . . . . . . 156 8.17 Contact pressure field at t = 2.5 ms for wet and dry smooth roads. . . . . . . 157 8.18 Time history of the wheel center vertical displacement for wet and dry smooth roads. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157 8.19 Aquaplaning configuration with a single APR box. . . . . . . . . . . . . . . . 158 8.20 Time history of the vertical fluid force with and without APR, smooth road. 159 8.21 Time history of wheel center vertical displacement with and without APR, smooth road. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 159 8.22 Time history of the vertical solid force on tire with and without APR, smooth Velocity field at different instants with (left) and without (right) APR, smooth road. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.1 Time history of the ratio between E Interf ace and E T OT for different stabilization coefficients, force conservative formulation. . . . . . . . . . . . . . . . . A.2 Time history of the ratio between E Interf ace and E V for different stabilization coefficients, force conservative formulation. . . . . . . . . . . . . . . . . . . . A.3 Pressure field comparison for different loading formulations (columns) at different instants (lines), HHT scheme with α s =-0.05. . . . . . . . . . . . . . A.4 Time history of the ratio between E Interf ace and E T OT for different stabilization coefficients, energy conservative formulation. . . . . . . . . . . . . . . . A.5 Time history of the ratio between E Interf ace and E V for different stabilization coefficients, energy conservative formulation. . . . . . . . . . . . . . . . . . . A.6 Time history of the ratio between E Interf ace and E V for different boundary conditions and fluid loading, HHT scheme with α s = -0.1. . . . . . . . . . . A.7 Time history of the atio between E Interf ace and E T OT for different boundary conditions and fluid loading, HHT scheme with α s = -0.1. . . . . . . . . . . A.8 Time history of the ratio between E Interf ace and E T OT for different stabiliza-Physical and numerical parameters of the hydrostatic water column in the elastic plate test case. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Physical and numerical parameters. . . . . . . . . . . . . . . . . . . . . . . . 7.3 Physical and numerical parameters. . . . . . . . . . . . . . . . . . . . . . . . 8.1 Aquaplaning simulation parameters. . . . . . . . . . . . . . . . . . . . . . . . 8.2 Averaged forces and solid-solid contact surfaces, smooth road. . . . . . . . . 8.3 Aquaplaning simulation parameters. . . . . . . . . . . . . . . . . . . . . . . .

	LIST OF TABLES
	List of Tables
	4.1

road. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160 8.23 Time history of the solid-solid contact surface with and without APR, smooth road. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160 8.24 Tire profiles at the left (top), central (middle) and right (bottom) rib crosssections at t = 2.5 ms with and without APR, wet smooth road. . . . . . . . 161 LIST OF FIGURES 8.25 tion coefficients. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.9 Time history of horizontal (top) and vertical (bottom) displacements of the gate extremity for different loading formulations, HHT scheme with α s = -0.1.172 A.10 Time history of horizontal force on the elastic gate for different loading formulations, Newmark scheme with α s = -0.3. . . . . . . . . . . . . . . . . . . A.11 Pressure field comparison for different loading formulations (columns) at different instants (lines), HHT scheme with α s =-0.05. . . . . . . . . . . . . . B.1 Pressure field comparison for different time step ratios p at different instants, CSS procedure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.2 Pressure field comparison for different time step ratios p at different instants, CSS procedure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.3 Pressure field comparison for different time step ratios p at different instants, CPS procedure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1 Physical and numerical parameters for the deformable beam impact. . . . . . 5.2 Physical and numerical parameters of the Antoci dam break. . . . . . . . . . 7.1 Physical and numerical parameters for the 3D dam-break through an elastic gate. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2

Table 5 .

 5 5.16). 2: Physical and numerical parameters of the Antoci dam break.

	L	79 mm
	H e l Poisson coefficient ν ρ rubber ρ water	0.14 m 5 mm 0.1 m ≈ 0.5 1100 kg.m -3 1000 kg.m -3	∆X F E ∆Y F E R/∆X SP H C 0 CFL number	1.25 mm 1.975 mm 4.0 30.0 ms -1 0.375
	g	9.81m.s -2

Table 7 .

 7 2: Physical and numerical parameters.

	∆x F E	0.5 mm
	∆y F E	1.435 mm
	∆z F E	4.15 mm
	∆x SP H	1 mm
	R/∆x SP H	2.11
	C 0	15.0 m.s -1
	CFL	0.375
	number

Table 7 .

 7 3: Physical and numerical parameters.

		5 mm
	∆y F E	0.7175 mm
	∆z F E	4.15 mm
	∆x SP H	1 mm
	R/∆x SP H	2.11
	C 0	8.0 m.s -1
	CFL	0.375
	number

Table 8 .

 8 2: Averaged forces and solid-solid contact surfaces, smooth road.

			Averaged
	Ground vertical force contact surface
		(N)	(mm 2 )
	Dry	4703	9393
	Wet	4295 (-8.7%) 7899 (-15.9%)

Remerciements

Discussion

The aim of this chapter was to present the SPH method ingredients (from the basics to advanced adaptivity techniques) used at the LHEEA laboratory. The advantages and limitations of the method have also been discussed. It allows to then choose the most adaptable SPH scheme for fluid-structure interactions problems. The regularity of the pressure fields appear to be essential for the quality of the coupling strategy, which preferably leads to the choice of a Riemann-SPH scheme. For the following simulations, the choice of the boundary condition method or adaptivity techniques will depend on the problem to model. CHAPTER 5. ENERGY CONSIDERATIONS ON THE SPH-FE COUPLING STRATEGY

Energy balance in the δ-SPH scheme

The energy balance for two SPH schemes is provided in this section, with emphasis on irreversible internal and boundary energies. The contributions of each energy term will be expressed independently of each other. We start with the δ-SPH method, before its analogy with a Riemann-SPH scheme which is then used throughout this study. This section proposes a methodology to estimate all energy components using the discrete equations. The total energy can be expressed as the sum of three terms: potential, kinetic and internal energies.

E T OT = E P + E K + E I .

(5.4)

As proposed by [START_REF] Antuono | Energy balance in the δ-SPH scheme[END_REF] (see Section 2.3.2), the set of discrete equations is defined as:

)

(5.8)

Potential energy

The variation of potential energy is classically defined as:

where N F P is the total number of particles in the fluid domain.

Kinetic energy

The variation of kinetic energy per unit mass E m K is expressed as:

(5.10)

The right hand side terms in Eq. ( 5.10) correspond respectively to the power of volumetric forces, the power of pressure forces -V . gradP and the power of viscous constraints. In this study an inviscid fluid is considered (τ = 0). The variation of kinetic energy per unit mass therefore writes:

(5.11)

Summary

The energy components of the δ-SPH scheme have been presented separately in detail, highlighting the presence of a boundary term which represents the power yielded or received through the solid boundaries. Different dissipative contributions which take part in increasing the irreversible internal energy have also been outlined. This dissipative energy is critical for the stability of SPH simulations.

5.3

Energy balance in the Riemann-SPH scheme

Riemann-SPH scheme with local mass conservation

The simulations discussed in the next sections are performed using a Riemann-SPH scheme.

Here, the formulation proposed by LEDUC et al. ( 2010) is used, ensuring the local mass conservation. We show here that this scheme is very close to the δ-SPH one. The Riemann scheme considered is defined by the following system:

)

dm i dt = 0, (5.23)

2P e ∇ i W ij ω j ω i , (5.24)

where P e and v e are the linearized Riemann problem solutions. Eq. ( 5.23) leads to:

(5.26)

These mass and momentum conservation equations are very similar to the ones in δ-SPH scheme (Eq. (5.6) and Eq. (5.7). Using:

(5.27) Π ij = (P i + P j -2P e ), (5.28) the mass and momentum conservation equations write:

CHAPTER 5. ENERGY CONSIDERATIONS ON THE SPH-FE COUPLING STRATEGY

(5.30)

Eq. (5.19) is also verified here. The analogy with Eq. ( 5.18) and Eq. ( 5.20) for internal irreversible and boundary energies is proposed as: (5.31)

(5.32)

These formulations are only valid for the ghost particle method. Using the Normal Flux Method presented in Section 2.5.1.2, energy equations need some rearrangements.

Normal Flux Method

Kinetic energy

The variation of the kinetic energy is still obtained by multiplying the momentum equation by the velocity of the particle i. Then a summation is done over all the fluid particles:

(5.33)

Internal energy

Following the same process:

Once more the term of Eq. ( 5.34) can be splitted into in two terms:

-

. (5.35) Finally: (5.36) Instabilities which appeared at low dissipation levels in the solid with the averaging pressure technique (Fig. 5.11 and Fig. 5.12) are absent here. The coupling is stable for any α s parameter values and whatever the time scheme used (HHT or Newmark). High frequency oscillations vanished from the fluid force time history for low dissipation values in the HHT scheme (see Fig. 5.29, where the results using the averaging pressure is recalled). These schemes are easily implementable and do not require any other change. Note that it is necessary to preserve a solid time step in agreement with the physics (i.e. period of natural vibrations). Concerning the computational costs, the speedup S p of these new algorithms can be defined as: S p = t original scheme t improved by p ratio .

Ghost particle method NFM

(6.6)

For each procedure the speedup can be expressed using t F E and t SP H , being the CPU STRATEGY Chapter 7

3D validation of the SPH-FE coupling strategy

All previous simulations were performed on 2D test cases. However, it seems important to investigate the relevance of our SPH-FE model on 3D validation test cases before considering more complex problems such as the 3D aquaplaning simulation. To our knowledge, the only validation case involving a 3D SPH-FE coupling has been proposed by [START_REF] Gilbert | Accelerating an SPH-FEM Solver using Heterogeneous Computing for use in Fluid-Structure Interaction Problems[END_REF].

In this chapter, we propose to validate our FSI strategy on three different experimental test cases. The first one corresponds to the 3D extension of the previous dam-break flow. Then two FSI problems considering confined sloshing flows with elastic beams are investigated. Note that the simulations are performed here with the Abaqus software from Simulia instead of Code Aster. Furthermore, the HHT scheme with α s = -0.3 is systematically used since the Newmark scheme is not available in Abaqus.

Dam-break flow through an elastic gate

The first validation of the 3D SPH-FE coupling strategy applies to the dam-break flow through an elastic gate (see Fig. Here, fluid particles are distributed uniformly to represent the fluid. The FE mesh is composed of 4 width-wise, 4 depth-wise and 40 length-wise elements. The rubber behaviour is still considered as incompressible with a Poisson coefficient close to 0.5. The non-linear strain-stress curve presented in Fig. 5.16 is still used to model the rubber gate behaviour. The boundary conditions are modelled with ghost particles. Furthermore, the CPS procedure with p = 2 is retained as the best compromise among accuracy, stability and computational time. A ratio p = 8 was retained in the previous chapter but using a finer fluid spatial resolution, ∆x SP H = 0.25 mm instead of ∆x SP H = 1 mm here. The motion is imposed to the wheel center through horizontal and rotation velocities, corresponding to 80 km.h -1 . Note that no friction is considered in all the following simulations. Furthermore, a HHT scheme is systematically used for the solid implicit time integration, with a damping value α s = -0.3.

The initialization of the solid part is performed through the following four steps:

1. Inflation at 0.2 MPa, static.

2. Vertical position setting in a first approximation with respect to the ground, static.

3. Vertical position setting by applying a load of 337 daN down to the wheel center, corresponding to the car weight on a single wheel, static.

4. Horizontal translation and rotation imposed to the wheel center. Note that no vehicle damper is considered here. The wheel is free to move in the vertical direction, dynamic.

Note that the contact between the tire and the road is handled by a contact pair algorithm from Abaqus software using a penalty method. In this technique the contact force is proportional to the tire penetration into the ground according to some degrees of interpenetration. We refer the reader to the documentation of Abaqus software SIMULIA (2012) for more informations about this contact management.

Numerical results

The global characteristics of the aquaplaning simulations are detailed in the previous section. Here, preliminary numerical results are presented. Fig. 8.3 shows the complexity of the free surface at t = 10 ms (the tire interface is not represented), outlining the complexity of the free surface deformations.

Appendix A

Numerical results for different fluid loading formulations using the NFM method A.1 The deformable beam impact Aquaplaning -------------------------------------------------------------------------------------------------------------------------------------

Abstract

The aquaplaning problem has been the topic of simulation works emphasizing its complexity: fluidstructure interactions, structures modelling, materials involved, contact with asphalt and the complexity of the resulting fluid flow (extremely complex interface, drying up the road, ventilation, possible development of turbulence and cavitation). As additional difficulty, the tire is a highly deformable body and fluid-structure interaction effects should be considered, leading to a challenging problem for the numerical modelling. Then Michelin, Ecole Centrale Nantes and NextFlow Software have recently tested the ability of the SPH solver developed by the two latter to classify tires based on their surface structure geometries, without considering the gas phase.

In this context, the interest of SPH for modelling efficiently the aquaplaning flow has been underlined. The meshless and Lagrangian feature of SPH naturally avoid the problem of fluid/solid grid compatibility. The other significant advantage of the SPH method, in this context, appears in its ability to be relatively easily coupled to with conventional Finite Element solvers. The aim of this work is threefold. First, qualify the SPH-FE coupling strategy, especially in terms of energy and then develop schemes to ensure a good compromise among stability, accuracy and computation time. Second, quantify the influence of different involved physical phenomena to determine which should be modelled. Finally, adapt SPH models to simultaneously consider different phenomena and to performe simulations of the complete problem.