Je Tiens

MASSIN, Directeur M Patrick

Stéphane Aubert

Je Voudrais

M Michel Lance

M Julien Leduc

Ingénieur M Hervé Ozdoba

Annabelle Zhang

Benoit Hennebelle

Guanyu Zhou

Feng Gao

Liuqing Wang

Gang Huang

Chang- Wei Zhou

Xu Chen

Guillaume Despres

Magdalena Neuhauser

Gilles-Alexis Renaut

Keywords: SPH method, finite element method, fluid-structure interaction, coupling algorithm, energy conservation at interface, numerical stability

Je tiens à adresser en premier lieu mes plus chaleureux remerciements à M. Alain COMBESCURE, Professeur de l'INSA de Lyon, le co-directeur de cette thèse, pour la confiance qu'il m'accordée en acceptant d'encadrer ce travail doctoral, pour ses précieux conseils scientifiques et analytiques, ainsi que pour sa passion de la recherche qu'il ma transmise pleinement. Tout au long de ces trois années, j'ai pu apprécier non seulement sa dimension scientifique, mais aussi sa non moins importante dimension humaine.

Je tiens à remercier du plus profond du coeur M. Francis LEBOEUF, Professeur de l'Ecole Centrale de Lyon, l'ancien directeur de ma thèse. Son décès, qui a eu lieu quelques mois avant la soutenance, m'a provoqué un choc émotif fort. Tout le monde a traversé un moment très difficile. Je n'oublierai jamais son enthousiasme et attachement pour la recherche scientifique, ses encouragements et ses conseils qui m'ont poussé énormément, son chapeau marron, son visage aimable, son élégance ... J'espère que ce manuscrit de thèse pourra être un bon cadeau pour Francis, et il sera un souvenir précieux pour moi.

Résumé

L'Interaction Fluide-Structure (IFS) est un sujet d'intérêt dans beaucoup de problèmes pratiques aussi bien pour les recherches académiques ainsi que pour les applications industrielles. Différents types d'approches de simulation numérique peuvent être utilisés pour étudier les problèmes d'IFS afin d'obtenir de meilleurs conceptions et d'éviter des incidents indésirables.

Dans ce travail, le domaine du fluide est simulé par une méthode hybride sans maillage (SPH-ALE), et la structure est discrétisée par la méthode d' Éléments Finis (EF). Considérant le fluide comme un ensemble de particules, on peut suivre l'interface entre le fluide et la structure d'une manière naturelle. Une stratégie de couplage conservant l'énergie est proposée pour les problèmes d'IFS transitoires où différents intégrateurs temporels sont utilisés pour chaque sous-domaine: 2 nd ordre schéma de Runge-Kutta pour le fluide et schéma de Newmark pour le solide. En imposant la continuité de la vitesse normale à l'interface, la méthode proposée peut assurer qu'il n'y a ni injection d'énergie ni dissipation d'énergie à l'interface. L'énergie de l'interface est donc nulle (aux erreurs de troncature près) durant toute la période de simulation numérique. Cette méthode de couplage assure donc que la simulation de couplage est numériquement stable en temps. Les expérimentations numériques montrent que le calcul converge en temps avec l'ordre de convergence minimal des schémas utilisés dans chaque sous-domaine.

Cette méthode proposée est d'abord appliquée à un problème de piston mono-dimensionnel. On vérifie sur ce cas qu'elle ne dégrade pas l'ordre de précision en temps des schémas utilisés. On effectue ensuite les études des phénomènes de propagation d'ondes de choc au travers de l'interface fluide-structure. Un exellent accord avec la solution analytique est observé dans les cas de teste de propagation d'onde en 1-D. Finalement, les exemples multi-dimensionnels sont présentés. Ses résultats sont comparés avec ceux obtenus par d'autres méthodes de couplage.

Mots-clés: méthode SPH; méthode d'éléments finis; intéraction fluidestrucutre; algorithme de couplage; conservation de l'énergie à l'interface; stabilité numérique. 6.9 Comparison between the analytical solution and the numerical results using the 1 st initial configuration with the implicit Newmark scheme. . . 6.10 Comparison between the analytical solution and the numerical results using the 1 st initial configuration with the explicit Newmark scheme. . . 6.11 Comparison between the numerical result and the analytical solution Let's consider the vibration of aircraft wings or turbomachinery blades, the flow of blood through arteries, the response of bridges and tall buildings to winds, etc. In all these phenomena a common point can be found out, which is the interaction between two different types of mediums, the fluid flows and the solid structures. The abbreviation FSI, which means Fluid-Structure Interaction, is widely used to denote this kind of phenomenon. In some cases, the interaction is only dedicated to the exchange of momentum, but not to the heat transfer or any other types. In such FSI problems, the movement and deformation of the structure are influenced by the fluid flow state, at the same time, the fluid flow is also affected by the position and movement of the solid structure. Concisely speaking, the interaction between the two mediums is mutual and simultaneous.

To understand the FSI phenomena, people try to construct some mathematical models and obtain the results by resolving analytically the appropriate differential equations. However, in most cases this method is quite difficult or even impossible to realize because of the complexity and nonlinearity of the differential equations. Sometimes it's hard to find the solutions, most of the time there exists no analytical solutions. In this case, numerical methods, such as Finite Element (FE) method or Finite Difference method (FD), are used to resolve numerically the complex and nonlinear differential equations.

Even though this makes the resolution possible, the numerical methods often suffer from other difficulties, for example, the numerical instability issues, the compatibility of different numerical methods used in each sub-domain, and the requirement of large computer resources etc. Hence, some approximative hypothesis have been made to simplify the resolution of FSI problems in realistic engineerings, such as neglecting the effect from one medium to the other or supposing that the interaction effect is static. Unfortunately, such simplification in conception may cause accidents or disasters, for example, the collapse of Tacoma Narrows Bridge due to the vibration induced by the surrounding air flows, and the fracture of water pipe caused by the propagation of shock waves (Figure 1.1), which is also called the Water Hammer problems in pipes [START_REF] Bergeron | Water hammer in hydraulics and wave surges in electricity[END_REF]56]. Hence, developing a robust numerical coupling method for FSI problems becomes an important issue for scientific research as well as realistic industrial engineering.

Overview of existing methods for FSI 1.2.1 Lagrangian vs. Eulerian points of view

Typically, in continuum mechanics, a physical problem can be described in two different ways: the Lagrangian approach and the Eulerian approach. In the Lagrangian approach, the physical variables, e.g. the density or the velocity etc., are defined at the infinitesimal material domain which is moving with the material point. Hence, this La-grangian description permits an easy tracking of free surfaces and interfaces between different mediums. In addition, the constitutive equations are always evaluated at the same material points, which is sensible for history-dependent materials. But it is hard to follow large distortions of the computational domain without remeshing operations. In the Eulerian approach, the variables are defined at the point which doesn't move in space as time varies, so the continuum moves with respect to the mesh. Large distortions in the continuum can be easily handled, but the drawback of this approach is the imprecise treatment of moving boundaries.

Based on the different characteristics of movement and state of solid and fluid mediums, it's often preferable to use Lagrangian approach in solid mechanics (the deformation apprears often moderate), and to use Eulerian approach in fluid mechanics.

However, a difficulty may occur for the treatment of mesh for the fluid sub-domain, as presented previously. Because, in Eulerian approach the fluid mesh is fixed in space, however, as the solid mesh is time-varying, there will be mesh incompatibility issues at the fluid-structure interface.

Various methods are used to overcome this difficulty, such as the Level Set (LS) method [START_REF] Osher | Fronts propagating with curvature dependent speed: algorithms based on Hamilton-Jacobi formulations[END_REF], the Volume of Fluid (VOF) method [START_REF] Hirt | Volume of fluid (VOF) method for the dynamics of free boundaries[END_REF] and the Immersed Boundary (IB) method [START_REF] Peskin | The immersed boundary method[END_REF] etc. These methods mentioned above are all based on the boundary tracking techniques on a fixed background mesh with the Eulerian approach for the fluid sub-domain. Nevertheless the most widely used method is named as the Arbitrary-Lagrangian-Eulerian (ALE) method [START_REF] Donea | Arbitrary Lagrangian-Eulerian Methods[END_REF], which combines the advantages of the above classical Lagrangian and Eulerian descriptions, while minimizing their respective drawbacks as far as possible.

In the ALE method, the nodes of the computational mesh may move with the continuum (Lagrangian description), or be held fixed in space as time varies (Eulerian description), or, move with an arbitrary velocity. When using this method to the fluid sub-domain, it allows us to easily track the fluid-structure interface, and to follow the large distortions or movements in the fluid continuum. However it requires additional resolutions for updating the fluid mesh. This will need more computational resources and sometimes lead to numerical instability issues, if the ALE mesh movement is not properly calculated.

In this work, the fluid sub-domain is discretized by a meshless particle approach in ALE description, and a Lagrangian mesh is used for the solid sub-domain. Since both the fluid and solid sub-domains are described in Lagrangian ways, it is easier and more natural to track the fluid-structure interface, and also to handle the free surfaces flows problems [START_REF] De Leffe | SPH modeling of shallow-water coastal flows[END_REF][START_REF] Maruzewski | SPH highperformance computing simulations of rigid solids impacting the free-surface of water[END_REF].

Monolithic vs. Partitioned algorithms

Generally the coupling strategies or algorithms of the two different physical domains can be classified into two categories [START_REF] Farhat | Robust and provably second-order explicit-explicit and implicit-explicit staggered timeintegrators for highly non-linear compressible fluid-structure interaction problems[END_REF]: monolithic and partitioned procedures.

In monolithic procedures [START_REF] Blom | A monolithical fluid-structure interaction algorithm applied to the piston problem[END_REF][START_REF] Michler | A monolithic approach to fluid-structure interaction[END_REF], the fluid and solid governing equations are solved simultaneously or synchronously. Hence, in most cases, the resolution has to be carried out in one numerical simulation solver or software. Preservation of accuracy and numerical stability is the most important advantage of this sort of method. This monolithic method is however difficult to handle if the fluid and the structure are described by two very different solvers.

Partitioned procedures [26,[START_REF] Matthies | Partitioned strong coupling algorithms for fluid-structure interaction[END_REF][START_REF] Piperno | Partitioned procedures for the transient solution of coupled aeroelastic problems -part ii: Energy transfer analysis and threedimensional applications[END_REF][START_REF] Wall | A strong coupling partitioned approach for fluid-structure interaction with free surfaces[END_REF] can overcome this limit, since the equations of fluid and structure are integrated in time separately and alternatively, and different numerical methods can be used in each sub-domain. However, in a partitioned scheme, the interface conditions are often enforced asynchronously, which may lead to a time lag between the time integration of the fluid and the structure [START_REF] Blom | A monolithical fluid-structure interaction algorithm applied to the piston problem[END_REF]. As a consequence, this partitioned method is sometimes energy increasing [START_REF] Michler | A monolithic approach to fluid-structure interaction[END_REF]. Considering the time discretization (different time schemes, different time scales for the fluid and the structure), such approaches may suffer from a lack of accuracy, non-optimal convergence rate, and possible numerical instabilities.

In 2001, Piperno and Farhat [START_REF] Piperno | Partitioned procedures for the transient solution of coupled aeroelastic problems -part ii: Energy transfer analysis and threedimensional applications[END_REF] proposed a very important improvement of partitioned solutions: it allows to recover an optimal convergence rate with a control of the energy balance at the interface. More recently Farhat et al., [START_REF] Farhat | Time decomposed parallel time integrators: theory and feasibility for fluid structure and fluid structure applications[END_REF]22,[START_REF] Farhat | Provably second-order time-accurate loosely-coupled solution algorithms for transient nonlinear computational aeroelasticity[END_REF][START_REF] Farhat | Robust and provably second-order explicit-explicit and implicit-explicit staggered timeintegrators for highly non-linear compressible fluid-structure interaction problems[END_REF], have shown that if an accurate time integration of the moving grid is done and the position of the structure is predicted by a second-order time-integrator, then the resulting solution procedure is formally second-order accurate, even though numerical stability is not formally proved. In 2004, Michler, Hulshoff, van Brummelen and de Borst [START_REF] Michler | A monolithic approach to fluid-structure interaction[END_REF] compared partitioned and monolithic solution procedures for the simulation of fluid-structure interaction. They show that for the partitioned procedures, the time lag between the integration of the fluid and the structure implies that the interface conditions cannot be satisfied exactly. This induces an algorithmic energy production or dissipation at the interface, which may cause numerical instability and restrict the admissible time step size. Although stability and accuracy of partitioned schemes can be improved by means of predictions techniques, their error remains larger than for a monolithic solution method [START_REF] Michler | A monolithic approach to fluid-structure interaction[END_REF], for which the interaction of the fluid and the structure at the mutual interface is treated synchronously.

The approach presented in the work of this thesis is a monolithic solution procedure for the numerical simulation of Fluid-Structure Interactions as a fully coupled problem.

Motivations

The laboratory LMFA (Laboratoire de Mécanique des Fluides et d'Acoustique) has developed with ANDRITZ Hydro since 10 years a numerical method based on a mesh-less technique with the aim of simulating and studying the free surface flows in hydraulic turbines, especially for the Pelton turbine which is a widely used water impulse hydraulic turbine. A typical Pelton turbine's configuration is shown in Figure 1.2. The water jet impacts on the turbine blade (Figure 1.3) with a very high speed so that the turbine wheel can be driven to rotate. As a consequence, the turbine rotates with a high rotating speed under a periodic loading charge by the water jet, which may induce the vibration or the fracture of the turbine blade. The rupture of blade will lead to a decrease of the turbine efficiency, or more dangerous, the cracked blade may crash the other blades or the turbine case. Some one-way coupling procedures (fluid-to-solid) may be used to simulate the water impacting effect on the blade [START_REF] Marongiu | Low mach number numerical schemes for the SPH-ALE method[END_REF][START_REF] Marongiu | Free surface flows simulations in Pelton turbines using an hybrid SPH-ALE method[END_REF][START_REF] Parkinson | From classical to novel flow simulation methods applied to Pelton rehabilitation projects[END_REF], as the deformation of the turbine's blade is normally very small, in such cases that the interaction effect due to the structural deformation is negligible.

Besides, lots of investigations in contact mechanics and coupling methods have been performed in the laboratory LaMCoS (Laboratoire de Mécanique des Contacts et des Structures).

Based on the knowledge of FSI methods and on the improvement of the SPH method to compute the fluid pressure fields, this PhD work is aimed at proposing and validating a two-way coupling strategy without time lag for numerical simulation of the fluidstructure interaction problems with large deformation or displacement of the structure in free surface flows where the usage of the SPH method is relevant.

Numerically, to track the fluid-structure interface and handle the free surface flows, we choose to use the Smoothed Particle Hydrodynamics (SPH) method to discretize the fluid sub-domain. And the semi-discrete fluid equations will be integrated in time with the 2 nd order Runge-Kutta scheme. For the structure, the Finite Element (FE) method will be applied to simulate the solid sub-domain with a Lagrangian mesh. The Newmark scheme is used to realize the time integration of the solid semi-discrete equations. The Finite Element numerical simulation is carried out in an open-source FE software which is called Code Aster developed by EDF R&D.

In this context, we developed a coupling strategy for numerical simulation of Fluid-Structure Interaction (FSI) where different time and space schemes are used for the two sub-domains. This proposed coupling strategy obtains the interface status by resolving a system of equations at the fluid-structure interface, which is naturally compatible with the fluid and solid sub-domains. From the viewpoint of energy, this presented coupling approach can conserve the interface energy as time advances. Hence, this coupling method can preserve the numerical stability of the used numerical schemes in each subdomain without degrading the optimal order of accuracy in time.

Structure of the thesis

This thesis is organized as follows: Chapter 2 presents the physical problems and the mathematical description tools, and then the governing equations used for the fluid and structure sub-domains; the numerical time and space schemes will be presented in Chapter 3 and Chapter 4; Chapter 5 shows the interface condition to be imposed in order to conserve the energy of the coupled system, based on which we obtain the system of equations at the interface to solve; to prove the quality of the proposed coupling method, some numerical examples will be presented in Chapter 6; and finally, the conclusions and perspectives will be given in Chapter 7.

Chapter 2

Physical problems and governing equations

Generally, a physical problem can be described by a set of mathematical equations, which are often some partial differential equations. By solving these governing equations, if possible, with the appropriate Boundary Conditions (BC) and Initial Conditions (IC), one can evaluate analytically the governed physical problem. However, the reliability and accuracy of the results will depend on the mathematical model, i.e. the chosen governing equations.

This chapter begins with a brief presentation of the continuum mechanics, in which the physical variables are expressed and the governing equations are constructed. Then a theorem used to derive the conservation equations, the Reynolds Transport Theorem (RTT), will be presented secondly. The third part will describe the physical problems investigated in the fluid and structure sub-domains. The last part will give the corresponding governing equations for each sub-domain.

Continuum mechanics

Continuum mechanics is a branch of mechanics that deals with the analysis of the kinematics and mechanical behavior of materials modeled as a continuous domain, rather than of the interaction effects between the molecules, such as molecular dynamics etc. The objective of continuum mechanics is to provide models for the macroscopic be-havior of fluids, solids. In the following, some specialized terminology in continuum mechanics will be presented.

Initial and current configurations

Consider a body in an initial state at the time t = 0 s, as shown in Figure 2.1. The material domain of the initial state is called the initial configuration and denoted by Ω 0 which is bounded by a closed surface Γ 0 . In many cases, it is also called the reference configuration, because the variables and equations are usually referred to this configuration. If in the initial state the body is not deformed, Ω 0 also denotes the undeformed configuration. Hence, if not specified otherwise, in most cases the initial configuration, the reference configuration and the undeformed configuration are identical. After a certain time, the body moves to a new position and possesses a new domain, which is called the current configuration denoted by Ω which is bounded by Γ, as shown in Figure 2.1. For deformable body, it is also called the deformed configuration.

Lagrangian and Eulerian coordinates

In Figure 2.1, the vector X denotes the initial position of the material point p in the reference configuration Ω 0 . Notice that X is independent of time, i.e. it will keep the same value as time varies and is written as

X = Xe x + Y e y + Ze z (2.1)
where X, Y and Z represent the scalar components of the vector X in the three directions, and, e x , e y and e z are the unit vectors in the 3-dimensional space.

Because the vector variable X for a given material point does not change with time, it provides a label for each material point, and then the variables X are called the material coordinates or Lagrangian coordinates.

The current position of the material point p is written as

x = xe x + ye y + ze z (2.2)
where x denotes the current position vector as shown in Figure 2.1, with x, y and z being the scalar components of the vector x. The variables x give the spatial positions and are called the spatial coordinates or Eulerian coordinates.

An application can be used to describe the relation between X and x, which represents the motion of the body

x = Φ(X, t) (2.3)
or in the inverse direction

X = Φ -1 (x, t) (2.4)
where x is the position of the material point X at the time t.

Lagrangian and Eulerian descriptions

As presented previously in the part of introduction, a variable or function can be expressed with two different approaches: the Lagrangian description and the Eulerian description.

If a function is expressed as F (X, t), it depends on the material or Lagrangian coordinates and the time t, and this is called the Lagrangian description. The variable or function is defined at the initial coordinates which are not time-dependent. Whereas, if the function is expressed as f (x, t), it will depend on the spatial or Eulerian coordinates and the time t. This function is then defined at the current coordinates which are time-dependent, and it is called the Eulerian description.

One of the major differences between the two sorts of descriptions will appear when we carry out the time derivative of the function F (X, t) and f (x, t) with X held constant, which is called the material time derivative or total time derivative

dF (X, t) dt = ∂F (X, t) ∂t X + ∂F (X, t) ∂X t • dX dt = ∂F (X, t) ∂t X (2.5)
where dX/dt = 0 since X is held constant, i.e. not time-dependent, and we have

df (x, t) dt = ∂f (x, t) ∂t x + ∂f (x, t) ∂x t • dx dt = ∂f (x, t) ∂t x + ∂f (x, t) ∂x t • ∂Φ(X, t) ∂t X + ∂Φ(X, t) ∂X t • dX dt = ∂f (x, t) ∂t x + ∂f (x, t) ∂x t • ∂Φ(X, t) ∂t X (2.6)
where Φ is an application describing the body motion by x = Φ(X, t), and the time derivative holding the spatial coordinates constant is called the spatial time derivative.

Arbitrary Lagrangian-Eulerian description

To avoid the drawbacks of each of the two classical descriptions, a more generalized description method is needed, which is called the Arbitrary Lagrangian-Eulerian (ALE) approach. In the ALE description of motion, neither the material configuration nor the spatial configuration is taken as the reference. The physical variables and equations are expressed in a third domain which is named as the referential configuration where the coordinates are denoted by χ. With another application, this reference coordinates χ is related with the material coordinates X

χ = Ψ(X, t) (2.7)
and the motion of the computational mesh is described by

x = Φ(χ, t) (2.8)
Consider a function or variable F(χ, t) expressed with reference coordinates χ, its the material time derivative writes

dF(χ, t) dt = ∂F(χ, t) ∂t χ + ∂F(χ, t) ∂χ t • dχ dt = ∂F(χ, t) ∂t χ + ∂F(χ, t) ∂χ t • ∂Ψ(X, t) ∂t X (2.9)
Remark: Consider a physical quantity expressed in these three different description methods: f (x, t), F (X, t) and F(χ, t). Notice that f , F and F are three different functions. In the following of this thesis, the same symbol will be used to express the physical variable in different descriptions for the sake of brevity.

Displacement, velocity and acceleration

The displacement of a material point p is defined as the difference between the current position and the initial position, as shown in Figure 2.1

u(X, t) = x(X, t) -X = Φ(X, t) -X (2.10)
where u(X, t) denotes the displacement vector of the material point p at the time t.

The velocity v(X, t) is the rate of change in time of the position vector for a material point, i.e. the material time derivative of x(X, t) . With (2.5), we have

v(X, t) = dx(X, t) dt = d [u(X, t) + X] dt = ∂u(X, t) ∂t X ≡ u (2.11)
The acceleration a(X, t) is the rate of change of velocity of a material point, in the similar way, we have

a(X, t) = dv(X, t) dt = ∂v(X, t) ∂t X = ∂ 2 u(X, t) ∂t 2 X ≡ v ≡ ü (2.12)
Notice that the variables u(X, t), v(X, t) and a(X, t) are all expressed in Lagrangian description. In fact, they can also be expressed in Eulerian description, e.g. the velocity vector v(x, t). As shown in (2.6), the material time derivative of v(x, t) writes

dv(x, t) dt = ∂v(x, t) ∂t x + ∂v(x, t) ∂x t • dx dt = ∂v(x, t) ∂t x + ∂v(x, t) ∂x t • v(x, t) (2.13)
As presented in [START_REF] Belytschko | Nonlinear Finite Elements for Continua and Structures[END_REF], (2.13) can be also written as

dv(x, t) dt = ∂v(x, t) ∂t x + v(x, t) • ∇ x v(x, t) (2.14)
or shortly written as

dv dt = ∂v ∂t x + v • ∇ x v (2.15)
where ∇ x denotes the left gradient operator [START_REF] Malvern | Introduction to the mechanics of a continuous medium[END_REF] with respect to the spatial coordinates.

The left gradient of the velocity vector is a tensor which can be written in the matrix form

∇ x v =        ∂v x ∂x ∂v y ∂x ∂v z ∂x ∂v x ∂y ∂v y ∂y ∂v z ∂y ∂v x ∂z ∂v y ∂z ∂v z ∂z        (2.16)
where v x , v y and v z represent the three components of the vector v in the three directions.

Remark: More generally, a "scalar" is a zero order tensor, a "vector" is a first order tensor and, if not specified, a "tensor" is a second order tensor. The gradient of a "scalar" is a "vector", the gradient of a "vector" is a "tensor", and the gradient of a "tensor" is a third order tensor.

When the velocity is expressed with the reference coordinates χ, with (

+ ∇ x v(x, t) • ∂x ∂χ t • dχ dt = ∂v(χ, t) ∂t χ + ∇ x v(x, t) • ∂ Φ(χ, t) ∂χ t • ∂Ψ(χ, t) ∂t X (2.20)
Additionally, because x = Φ(X, t) = Φ(χ, t), the material velocity can be written as

v(x, t) = ∂Φ(X, t) ∂t X = ∂ Φ(χ, t) ∂t X = ∂ Φ(χ, t) ∂t χ + ∂ Φ(χ, t) ∂χ t • ∂Ψ(χ, t) ∂t X = v(χ, t) + ∂ Φ(χ, t) ∂χ t • ∂Ψ(χ, t) ∂t X (2.21)
where the time derivative of Φ(χ, t) holding χ constant represents the velocity of the computational mesh, and is denoted by v(χ, t). Hence, with (2.21), (2.20) becomes

dv(χ, t) dt = ∂v(χ, t) ∂t χ + ∇ x v(x, t) • [v(x, t) -v(χ, t)] (2.22)
or shortly written as

dv dt = ∂v ∂t χ + (v -v) • ∇ x v (2.23)
which is usually used to calculate the material time derivative of a variable expressed with the reference coordinates. It is very attractive since the motion of the mesh can be arbitrarily chosen. Hence, it allows to avoid the large distortion of the mesh, and to track the interface of different sorts of mediums. For any function f (χ, t), its material time derivative writes

df dt = ∂f ∂t χ + (v -v) • ∇ x f (2.

24)

Remark: The reference coordinates χ do not represent the position of the mesh points.

The mesh motion is the mapping from the referential configuration to the spatial configuration.

Reynolds Transport Theorem

To derive the basic conservation equations in continuum mechanics, a useful tool is often applied, which is called the Reynolds Transport Theorem (RTT), and is also known as the Leibniz-Reynolds transport theorem. [START_REF] Donea | Arbitrary Lagrangian-Eulerian Methods[END_REF] gives a generalized form of the RTT

∂ ∂t χ Ω f (x, t) dx = Ω ∂f (x, t) ∂t x dx + Γ [f (x, t) v(x, t)] • n dx (2.25)
where f (x, t) can be any function or quantity expressed in the spatial domain. Ω is an arbitrary time-varying domain bounded by the surface Γ, which is a function of the reference coordinates χ and the time t, hence we can write Ω = Ω(χ, t) and Γ = Γ(χ, t). Ω is also called the control volume, which moves with the velocity v. n is the outward-pointing normal vector at the boundary of Ω.

Here, χ can be arbitrarily chosen, which leads to the free choice of the velocity of the control volume. There exist two special cases, when χ = X, Ω = Ω(X, t) which represents a material domain, and hence the time derivative of the integral within this domain is time material derivative

∂ ∂t X Ω f (x, t) dx = Ω ∂f (x, t) ∂t x dx + Γ [f (x, t)v(x, t)] • n dx (2.26)
where v denotes the material velocity. When χ = x, then Ω = Ω(x, t) which is a spatial domain fixed in space, and

∂ ∂t x Ω f (x, t) dx = Ω ∂f (x, t) ∂t x dx (2.27)
where the surface integral term disappeared, since the velocity of the control volume is zero.

Governing equations for fluid and solid domains

The system of governing equations usually consists of the conservation laws of the mass, the momentum and the energy. All these conservation equations can be derived by using the RTT with χ being the material coordinates X, i.e. the control volume represents a material domain.

Conservation of mass Taking f (x, t) = ρ(x, t) in (2.26)

∂ ∂t X Ω ρ(x, t) dx = Ω ∂ρ(x, t) ∂t x dx + Γ [ρ(x, t)v(x, t)] • n dx (2.28)
where the integral of density ρ(x, t) over the material domain Ω equals the mass of this material domain. In this thesis, it is considered that within a material domain the mass can neither be created nor destroyed. Hence the material time derivative of the mass is zero, i.e.

Ω ∂ρ(x, t) ∂t

x dx + Γ [ρ(x, t)v(x, t)] • n dx = 0 (2.29)
with the Gauss's theorem, (2.29) can be rewritten as

Ω ∂ρ(x, t) ∂t x + ∇ x • [ρ(x, t)v(x, t)] dx = 0 (2.30)
since Ω can be any material domain, the differential conservation equation of mass is given as ∂ρ(x, t) ∂t

x + ∇ x • [ρ(x, t)v(x, t)] = 0 (2.31)
or, for the sake of brevity, (2.31) can also be written as

∂ρ ∂t x + ∇ x • (ρv) = 0 (2.32)
Conservation of linear momentum Taking f = ρv in (2.26)

∂ ∂t X Ω ρv dΩ = Ω ∂ρv ∂t x dΩ + Γ (ρv ⊗ v) • n dΓ (2.33)
where the left term represents the material time derivative of the linear momentum of the material domain Ω, which equals the sum of the forces exerted onto this material domain, i.e.

Ω ∂ρv ∂t x dΩ + Γ (ρv ⊗ v) • n dΓ = Ω ρb dΩ + Γ σ • n dΓ (2.34)
where b denotes the body force vector, σ the tensor of stress. Using once more the Gauss's theorem, (2.34) becomes

Ω ∂ρv ∂t x dΩ + Ω ∇ x • (ρv ⊗ v) dΩ = Ω ρb dΩ + Ω ∇ x • σ dΩ (2.35)
As the material domain Ω can be arbitrarily chosen, the differential conservation equation of linear momentum is given as

∂ρv ∂t x + ∇ x • (ρv ⊗ v) = ρb + ∇ x • σ (2.36) Conservation of energy Taking f = ρe + 1 2 ρv • v in (2.26) ∂ ∂t X Ω ρe + 1 2 ρv • v dΩ = Ω ∂ ∂t x ρe + 1 2 ρv • v dΩ + Γ ρe + 1 2 ρv • v v • n dΓ (2.37)
where e denotes the internal energy per unit volume. The material time derivative of the total energy over the material domain Ω equals the power of the volume and surface energy sources. Hence the balance of the energy can be expressed by

Ω ∂ ∂t x ρe + 1 2 ρv • v dΩ + Γ ρe + 1 2 ρv • v v • n dΓ = Ω v • ρb dΩ + Γ v • t dΓ + Ω ρs dΩ - Γ q • n dΓ (2.38)
where t denotes the external force vector exerted on the boundary of the material domain, ρs the heat source per unit volume, and q the heat flux per unit area.

As mentioned previously, the interaction between the two physical domains consists only in the exchange of the linear momentum, hence the thermal phenomena are neglected, i.e. there is no heat transfer at the interface, or heat creation inside the material bodies. As a consequence, in this thesis the equations of energy will not be used.

The fluid sub-domain

In fluid mechanics, the partial differential equations (2.32) and (2.36) are easily to be implemented in the framework of Finite Difference (FD) method. For the reason of the used numerical method for the fluid sub-domain, a integral form of the conservation equations is given by [START_REF] Donea | Arbitrary Lagrangian-Eulerian Methods[END_REF].

Consider that at the instant t, the control volume Ω and the material domain Ω occupy the same spatial volume, hence Ω ∂ρ(x, t) ∂t

x dx = Ω ∂ρ(x, t) ∂t x dx (2.39)
then with (2.25) and (2.26)

∂ ∂t X Ω f dΩ = ∂ ∂t χ Ω f dΩ + Γ≡ Γ [(v -v)f ] • n dΓ = ∂ ∂t χ Ω f dΩ + Ω ∇ x • [(v -v)f ] dΩ (2.40)
which represents an integral form of the ALE differential equation (2.24). Notice that in (2.40) the material time derivative (the left term) is identical with the one in (2.26).

Hence, when taking f = ρ f and f = ρ f v f , we have The water is chosen as the investigated medium for the study of the FSI effects, which is considered in this work as a non-viscous quasi-incompressible fluid. Because the viscous effect is neglected, the stress tensor reduces to

         ∂ ∂t χ Ω ρ f dΩ + Ω ∇ x • ρ f (v f -v) dΩ = 0 ∂ ∂t χ Ω ρ f v f dΩ + Ω ∇ x • ρ f v f ⊗ (v f -v) dΩ = Ω ρ f b dΩ + Γ σ f • n dΓ
σ f = -p f I (2.42)
where p f is the fluid pressure and I denotes the unit tensor. With the Gauss's theorem, (2.41) becomes

         ∂ ∂t χ Ω ρ f dΩ + Ω ∇ x • ρ f (v f -v) dΩ = 0 ∂ ∂t χ Ω ρ f v f dΩ + Ω ∇ x • ρ f v f ⊗ (v f -v) + p f I dΩ = Ω ρ f b dΩ (2.43)
In addition, the fluid flows are considered as isothermal-barotropic, which means that the fluid temperature is constant and the fluid pressure is only a function of the fluid density (or vice versa). In 1966 [START_REF] Macdonald | Some simple isothermal equations of state[END_REF], J.R. Macdonald proposed some simple isothermal equations of state. In this thesis, the modified Tait equation is chosen as the equation of state for the water

p f = B   ρ f ρ ref f ζ -1   (2.44)
where The equations (2.43) and (2.44) constitute the governing equations for the fluid subdomain, which are to be discretized in space and integrated in time.

B = ρ ref f (c ref f ) 2 /ζ

The solid sub-domain

Equations (2.32) and (2.36) can also be used for the solid sub-domain as the governing equations, where the physical variables are expressed in Eulerian description.

Eulerian description The conservation equations (2.32) and (2.36) are written in the conservative form, which is often used in fluid mechanics to verify more accurately the conservation properties. In solid mechanics, another form is preferred

     dρ s dt + ρ s ∇ x • v s = 0 ρ s dv s dt = ρ s b + ∇ x • σ s (2.45)
which is mathematically identical with (2.32) and (2.36). The total time derivative formulation is given in (2.15). If the solid is incompressible, the incompressibility condition is given by

∇ x • v s = 0 (2.46)
which is also applicable for the fluid.

Remark: An incompressible material or medium signifies that the material or total time derivative of the density is zero, but not the local or spatial time derivative.

This system of governing equations in Eulerian description is usually used with the Finite Element (FE) method in Updated Lagrangian Formulation (ULF). One important advantage of this ULF is that the stress and strain are defined in the spatial domain, which means that they both have the intuitive physical senses. Hence the constitutive models are easily to be implied in the numerical simulation.

As presented previously, a physical variable can also be expressed with the Lagrangian description approach, i.e. with the material coordinates X.

Lagrangian description Here the governing equations in Lagrangian description way are given directly. More details can be found in [START_REF] Belytschko | Nonlinear Finite Elements for Continua and Structures[END_REF].

     ρ s (X, t)J(X, t) = ρ 0 s (X) ρ 0 s (X) ∂v s (X, t) ∂t X = ∇ X • P + ρ 0 s (X)b (2.47)
where ρ 0 s (X) denotes the initial solid density which depends on only the material coordinates X, J(X, t) the determinant of the Jacobian matrix of the motion Φ(X, t), ∇ X the left gradient operator with respect to the material coordinates, and P is the nominal stress tensor defined in the initial configuration Ω 0 .

The corresponding formulation of FE method is called the Total Lagrangian Formulation (TLF), in which the weak form of (2.47) is obtained by integration over the initial configuration.

More details about the two formulations ULF and TLF will be offered in Chapter 4 which presents the numerical method for the solid sub-domain.

Remark: These two formulations ULF and TLF can be transformed from one to the other, i.e. ULF and TLF are theoretically identical. Besides, both of the two formulations use the Lagrangian meshes.

Chapter 3 Smoothed Particle Hydrodynamics numerical method for fluid

Various numerical methods can be used to simulate numerically the movement of the fluid, e.g. the Finite Volume (FV) method, the Finite Difference (FD) method [START_REF] Hirsch | Numerical computation of internal and external flows[END_REF][START_REF] Hirsch | Numerical computation of internal and external flows[END_REF] etc. Because of the difficulties of handling the motion of the computational mesh to track the fluid-structure interface, particle methods are a solution in the case of free surface flows with moving solid boundaries.

In this chapter, a well known particle method, the Smoothed Particle Hydrodynamics (SPH) method will be presented first, followed by a special formulation developed by J.P. Villa [START_REF] Vila | On particle weighted methods and smooth particle hydrodynamics[END_REF]. The treatment of the solid wall boundary condition will be presented in the second part. The third part synthesizes the semi-discrete equations for the fluid sub-domain. Finally, the Runge-Kutta scheme is shown as the used time integrator.

Smoothed Particle Hydrodynamics method

The Smoothed Partical Hydrodynamics (SPH) method was originally proposed in 1977 by R.A. Gingold and J.J. Monaghan [START_REF] Gingold | Smoothed particle hydrodynamics: Theory and application to non-spherical stars[END_REF], and independently by L.B. Lucy [START_REF] Lucy | Numerical approach to testing the fission hypothesis[END_REF] for modeling astrophysical phenomena, and later widely extended for applications to continuum solid and fluid mechanics.

In the SPH method, the physical domain is discretized into a set of material particles which move with respect to each other. As a Lagrangian mesh-less approach, the SPH method possesses several features, such as being natural to track the interface of different materials, and easy to handle the large deformation of the investigated medium. When simulating the free surface flows, there is no need to put the particles in the air, because the zero pressure condition at the free surface is automatically imposed. This reduces the time for numerical simulation. For these reasons, the SPH method is chosen as the numerical method to simulate the fluid sub-domain for the fluid-structure interaction.

The basis of the SPH method

Essentially the SPH method is a technique for approximating the governing equations in continuum mechanics. The basic idea of this method is to represent the value of a function of a certain particle i by using the values of other particles which are located sufficiently near to this particle i.

The concept of the SPH method starts with the integral representation of a function f (x)

f (x) = Ω f (x ′ )δ(x -x ′ ) dx ′ (3.1)
where x and x ′ denote two independent spatial coordinates in the space Ω = R d with d being the number of the dimensions, and δ(xx ′ ) is the Dirac delta function defined by

δ(x) = + ∞ , x = 0 0 , x = 0 (3.2)
which is also constrained to satisfy the identity condition

+∞ -∞ δ(x) dx = 1 (3.3)
In the SPH method, this delta function is replaced by a smoothing function W (xx ′ , h) which is called the the smoothing kernel function or kernel function. Hence, f (x) is approximated by

f (x) ≃< f (x) >= D f (x ′ )W (x -x ′ , h) dx ′ (3.4)
where h is the smoothing length which defines the influence area D of the smoothing kernel function W . In the SPH convention, < f (x) > denotes the kernel approximation of the function f (x).

As < f (x) > is only an approximation of f (x), the choice of the kernel function is very important. Some basic conditions for the kernel function W are presented in the following:

◮ Normalization condition D W (x -x ′ , h) dx ′ = 1 (3.5) ◮ Symmetry condition W (x -x ′ , h) = W (x + x ′ , h) ∇ x ′ W (x -x ′ , h) = -∇ x ′ W (x + x ′ , h) (3.6)
◮ Delta Dirac function property

lim h→0 W (x -x ′ , h) = δ(x -x ′ ) (3.7)
◮ Compact condition

W (x -x ′ , h) = 0 when |x -x ′ | > κh (3.8)
where κ is a constant defining the non-zero area of the kernel function, which is called the kernel ratio in this thesis.

An example of the kernel function is shown in Figure 3.1, in which i is the investigated fluid particle, j is one of the neighbor particles. The coordinates of the fluid particle i and j are x and x ′ , respectively, and r = xx ′ . Finally κh is the radius of the influence area D i of the kernel function, which is bounded by ∂D i .

If the spatial integral over the influence area (the support domain D i ) is approxi- mated by the sum carried out over all the neighbor particles, then (3.4) becomes

f (x i ) ≃< f (x i ) > = D f (x ′ )W (x i -x ′ , h) dx ′ ≃ j∈D i f (x j )W (x i -x j , h)ω j (3.9)
which is called the particle approximation, where ω j denotes the volume of fluid particle j in 3-dimensional cases. For the sake of brevity, (3.9) is often written as

f i ≃ j∈D i f j W ij ω j (3.10)
where

f i = f (x i ) and W ij = W (x i -x j , h).
These kernel and particle approximations can also be applied to the gradient or divergence of a function, i.e.

∇ x f (x) or ∇ x • f (x) when f (x) is a vector or tensor field. Similarly to (3.4), the kernel approximation of ∇ x f (x) is ∇ x f (x) ≃< ∇ x f (x) >= D ∇ x ′ f (x ′ )W (x -x ′ , h) dx ′ (3.11)
where ∇ x ′ denotes the gradient operator with respect to x ′ .

With the help of the Gauss's theory, we can obtain

∇ x f (x) ≃< ∇ x f (x) > = D ∇ x ′ f (x ′ )W (x -x ′ , h) dx ′ = ∂D [f (x ′ )W (x -x ′ , h)]n dx ′ - D f (x ′ )∇ x ′ W (x -x ′ , h) dx ′ (3.12)
where n is the out-pointing normal vector of the boundary ∂D. If the support domain D is not truncated by the solid wall or the free surface, the integral over ∂D will dismiss, since ∀x ′ ∈ ∂D, W (xx ′ , h) = 0. Hence, (3.12) becomes

∇ x f (x) ≃< ∇ x f (x) > = D ∇ x ′ f (x ′ )W (x -x ′ , h) dx ′ = - D f (x ′ )∇ x ′ W (x -x ′ , h) dx ′ (3.13)
And finally, with the particle approximation, the gradient of a function is approximated by

∇ x f (x)| x=x i ≃ < ∇ x f (x) >| x=x i ≃ - j∈D i f (x j )∇ j W (x i -x j , h)ω j = - j∈D i f j ∇ j W ij ω j (3.14)
where ∇ j is the gradient operator with respect to x j , and because of the symmetry property of the kernel function W , we have

∇ j W ij = -∇ i W ij . Then (3.14) becomes < ∇ x f (x) >| x=x i ≡ [∇ x f (x)] i ≃ j∈D i f j ∇ i W ij ω j (3.15)
Remark: The interesting point of the SPH method is that the derivation is all carried out on the known kernel function. (3.15) is for the fluid particle whose support domain is not truncated by the solid boundary or the free surface. The truncation problem will be presented in Section 3.2.

Consistency of the SPH method

Because of the kernel and particle approximation, generally the SPH method using (3.10) and (3.15) is not consistent even of zero order, if no correction techniques are applied.

A numerical method is said to be p th order consistent means that it can represent exactly a p th order polynomial field [START_REF] Belytschko | Meshless methods: An overview and recent developments[END_REF]. For example, if the numerical approximation of a gradient function is zero order consistent, then the gradient of a spatially constant function is supposed to be zero. However, the gradient calculated by (3.15) is generally non-zero, even if the kernel function is symmetric, because the fluid particle can be irregularly distributed in space. Similarly, for a linear function, the calculated gradient should be a constant which equals the rate of slope, if the used numerical method is first order consistent.

To make the SPH method consistent, several methods have been proposed, e.g. in [START_REF] Monaghan | An introduction to SPH[END_REF] Monaghan introduced a symmetrization procedure in which he assumed that

f (x i ) j∈D i ∇ j W (x i -x j , h)ω j ≡ 0 (3.16)
although this is only validate for a regular distribution of fluid particles which are far from any boundary. Adding (3.16) to (3.15) gives

[∇ x f (x)] i ≃ j∈D i (f j -f i )∇ i W ij ω j (3.17)
which can theoretically represent a gradient function, and numerically ensure the zero order consistency.

Based on (3.17), in 1996 Johnson et al. [START_REF] Johnson | SPH for high velocity impact computations[END_REF] and Randles et al. [START_REF] Randles | Smoothed Particle Hydrodynamics: Some recent improvements and applications[END_REF] introduced independently a normalization correction method to ensure a first consistency for the approximation of gradient. In this method, taking a linear function f (x) = b • x + a, a second order tensor B is used to correct the gradient operator for the particle i

[∇ x f (x)] i ≃ B i • j∈D i [b • (x j -x i )] ∇ i W ij ω j = B i • j∈D i ∇ i W ij ⊗ (x j -x i )ω j • b ≡ I • b (3.18)
hence, with

B i = j∈D i ∇ i W ij ⊗ (x j -x i )ω j -1 (3.19)
the new normalized gradient operator will possess a first order consistency.

Besides, other methods are also used to improve the accuracy of the SPH method, such as the Moving Least Squares (MLS) method [START_REF] Belytschko | Element-free Galerkin methods[END_REF]9,[START_REF] Dilts | Moving-least-squares-particle hydrodynamics I: Consistency and stability[END_REF][START_REF] Nayroles | Generalizing the finite element method: diffuse approximation and diffuse elements[END_REF], the Reproducing Kernel Particle (RKP) method [START_REF] Liu | Reproducing kernel particle methods[END_REF] etc. An analysis of truncation error for SPH method is shown in the work of N. J. Quinlan [START_REF] Quinlan | Truncation error in meshfree particle methods[END_REF].

SPH method in fluid mechanics

The most common SPH method in fluid mechanics uses Lagrangian description. However, the spatial gradient is carried out with respect to the Eulerian coordinates. The governing equations for the fluid domain are written in differential form

     dρ f dt = -ρ f ∇ x • v f ρ f dv f dt = ρ f b + ∇ x • σ f (3.20)
where σ f = -p f I for inviscid flows, and b is usually the gravity vector g.

Applying the SPH method (3.17

) to (3.20) gives [59]            dρ i dt = ρ i j∈D i (v j -v i ) • ∇ i W ij ω j ρ i dv i dt = ρ i g i - j∈D i (p j -p i )∇ i W ij ω j (3.21)
where the second equation is the momentum equation for the fluid particle i. As presented by Monaghan in [START_REF] Monaghan | Smoothed Particle Hydrodynamics[END_REF][START_REF] Monaghan | Smoothed particle hydrodynamics[END_REF], in this equation the estimated pressure gradient vanishes exactly when the pressure is constant in space, which is an advantage that leads to a zero order consistency. However, it does not conserve linear or angular momentum exactly, since the force on particle i owning to j is not equal and opposite to the force on j owning to i, which means

(p j -p i )∇ i W ij ω j ω i = -(p i -p j )∇ j W ij ω i ω j (3.22) since (p j -p i )∇ i W ij ω j ω i = (p i -p j )∇ j W ij ω i ω j , with ∇ i W ij = -∇ j W ij .
In this case, Monaghan proposed a symmetrizing procedure to the pressure gradient operator. The momentum equation rewrites

dv f dt = g - ∇ x p f ρ f (3.23) according to ∇ x p f ρ f = ∇ x p f ρ f + p f ρ 2 f ∇ x ρ f (3.24)
to which applying the initial SPH approximation (3.15), the momentum equation (3.23) becomes

dv i dt = g i - j∈D i m j p j ρ 2 j + p i ρ 2 i ∇ j W ij (3.25)
where m j = ρ j ω j denoting the mass of the fluid particle j. And now the conservation of momentum is ensured by

m i m j p j ρ 2 j + p i ρ 2 i ∇ j W ij = -m j m i p i ρ 2 i + p j ρ 2 j ∇ i W ij (3.26)
which describes the equilibrium of forces between the fluid particle i and j.

In addition, an artificial viscous term Π ij is often added into (3.25) to simulate shock phenomena, or simply to stabilize a numerical scheme

dv i dt = g i - j∈D i m j p j ρ 2 j + p i ρ 2 i + Π ij ∇ j W ij (3.27)
the expression of Π ij is given by Monaghan in [START_REF] Monaghan | Smoothed particle hydrodynamics[END_REF].

Finally, with the equation of motion for the fluid particle i, the semi-discrete equations for inviscid fluids are written as

                   dx i dt = v i dρ i dt = ρ i j∈D i (v j -v i ) • ∇ i W ij ω j dv i dt = g i - j∈D i m j p j ρ 2 j + p i ρ 2 i + Π ij ∇ j W ij (3.28)
Notice that (3.28) is often seen in the SPH literature, but it is not unique. Other schemes may also be used, if they satisfy basic needs as a numerical method. Some of them keep the conservation property, whereas some of them improve the accuracy.

Vila's formulation

In Chapter 2, the Reynolds Transport Theorem has been used to derive the governing equations (2.43) for the fluid sub-domain. This system of equations can be rewritten in a more concise way

∂ ∂t χ Ω Φ dΩ + Ω ∇ x • F dΩ = Ω S dΩ (3.29)
where Φ denotes the variables, F the flux tensor, and

F = ρ f (v f -v) when Φ = ρ f ρ f v f ⊗ (v f -v) + p f I when Φ = ρ f v f (3.30)
and S is the source term, when

Φ = ρ f , S = 0; when Φ = ρ f v f , S = ρ f g. Noting that Ω ∇ x • F dΩ = ω i [∇ x • F] i (3.31)
where [∇ x • F] i denotes the mean value of the divergence of F over the fluid particle volume ω i , which is then approximated by the SPH method (3.15)

[∇ x • F] i ≃ j∈D i F j • ∇ i W ij ω j (3.32)
With (3.31) and (3.32), (3.29) becomes

d(ω i Φ i ) dt + ω i j∈D i F j • ∇ i W ij ω j = S i ω i (3.33)
where Φ i denotes the mean value of Φ over the fluid particle i, and S i = ρ i g i . This semi-discrete equation (3.33) for the fluid sub-domain has the same form with the one obtained by J.P. Vila in [START_REF] Vila | On particle weighted methods and smooth particle hydrodynamics[END_REF]. As presented previously, this equation can not ensure the conservation property, because

ω i F j • ∇ i W ij ω j = -ω j F i • ∇ j W ij ω i (3.34) since usually F i = F j . Taking an unit function H(x) = 1, ∀x ∈ R d , hence ∇ x H(x) = 0. Then inserting F i • ∇ x H(x) into (3.29) gives ∂ ∂t χ Ω Φ dΩ + Ω (∇ x • F + F i • ∇ x H) dΩ = Ω S dΩ (3.35)
where F i can be seen as a constant tensor which equals the flux tensor for the fluid particle i.

Since F i • ∇ x H(x) = 0, (3.35) is theoretically identical with (3.29).
Using the SPH approximation to this additional term in (3.35)

Ω F i • ∇ x H dΩ = F i • Ω ∇ x H dΩ ≃ F i • ω i j∈D i H j ∇ i W ij ω j = ω i j∈D i F i • ∇ i W ij ω j (3.36)
where

H j = H(x j ) = 1.
With (3.36), the semi-discrete form for (3.35) writes

d(ω i Φ i ) dt + ω i j∈D i (F j + F i ) • ∇ i W ij ω j = S i ω i (3.37)
As mentioned by Vila in [START_REF] Vila | On particle weighted methods and smooth particle hydrodynamics[END_REF], (3.37) keeps exactly the conservation property, since

ω i (F j + F i ) • ∇ i W ij ω j = -ω j (F i + F j ) • ∇ j W ij ω i (3.38)
Applying (3.37) in one-dimensional cases, it turns out to be a unconditionally unstable numerical scheme [START_REF] Vila | On particle weighted methods and smooth particle hydrodynamics[END_REF]. Although it is possible to insert an artificial viscosity as in the Monaghan's formula (3.27), we choose not to use this method in this work.

Vila [START_REF] Vila | On particle weighted methods and smooth particle hydrodynamics[END_REF] proposed an alternative method based on (3.37), which makes use of the approximate Riemann solver to stabilize the numerical scheme. In this novel method,

(F i + F j ) is calculated by 2G E,ij = 2G E (Φ i , Φ j ) with G E,ij representing an intermedi-
ate flux tensor obtained by solving a Riemann problem between the two states of fluid Φ i and Φ j , so that (3.37) is written as

d(ω i Φ i ) dt + ω i j∈D i 2G E,ij • ∇ i W ij ω j = S i ω i (3.39)
with

G E,ij = ρ E,ij (v E,ij -v ij ) when Φ i = ρ i ρ E,ij v E,ij ⊗ (v E,ij -v ij ) + p E,ij I when Φ i = ρ i v i (3.40)
where ρ E,ij , p E,ij and v E,ij are the solutions of the Riemann problem and the arbitrary velocity

v ij = v(x ij , t) with x ij = (x i + x j )/2
. These solutions will be described in the next sub-section.

Riemann problem and Riemann solvers

The partial differential governing equations for inviscid fluid is given as

       ∂ρ ∂t x + ρ∇ x • v + v • ∇ x ρ = 0 ∂v ∂t x + v • ∇ x v + ∇ x p ρ = g (3.41)
which are identical with the conservative form (2.32) and (2.36) already presented in Chapter 1. Notice that, in these two equations, the fluid variables are all expressed with the time t and the spatial coordinates x. And | x means that the partial time derivative are carried out with respect to x. In the following, | x will be eliminated for the sake of brevity.

Taking an example in two-dimensional space, (3.41) writes

               ∂ρ ∂t + ρ( ∂v x ∂x + ∂v y ∂y ) + (v x ∂ρ ∂x + v y ∂ρ ∂y ) = 0 ∂v x ∂t + (v x ∂v x ∂x + v y ∂v x ∂y ) + 1 ρ ∂p ∂x = g x ∂v y ∂t + (v x ∂v y ∂x + v y ∂v y ∂y ) + 1 ρ ∂p ∂y = g y (3.42)
which are the classical two-dimensional Euler equations in non-conservative form for inviscid flows [START_REF] Hirsch | Numerical computation of internal and external flows[END_REF][START_REF] Hirsch | Numerical computation of internal and external flows[END_REF].

The following procedure has been realized and presented in the PhD thesis of J.C. Marongiu [START_REF] Marongiu | Méthode numérique lagrangienne pour la simulation d'écoulements à surface libre -Application aux turbines Pelton[END_REF] in 2007 and the one of J. Leduc [START_REF] Leduc | Étude physique et numérique de l'écoulement dans un système d'injection de turbine Pelton[END_REF] in 2010.

To calculate G E,ij in (3.39), we are going to study the interaction effect between the fluid particle i and j as shown in Figure 3.2, where n ij denotes the unit direction vector defined as

n ij = (x j -x i )/|x j -x i |.
For each pair of fluid particles, such as i and j, let's construct a new orthogonal coordinate system (e ξ , e η ) where e ξ is parallel to n ij . Expressing the fluid variables in Figure 3.2: The Riemann problem between the fluid particle i and the particle j. this new coordinate system, the system of equations (3.42) becomes

               ∂ρ ∂t + ρ( ∂v ξ ∂ξ + ∂v η ∂η ) + (v ξ ∂ρ ∂ξ + v η ∂ρ ∂η ) = 0 ∂v ξ ∂t + (v ξ ∂v ξ ∂ξ + v η ∂v ξ ∂η ) + 1 ρ ∂p ∂ξ = g ξ ∂v η ∂t + (v ξ ∂v η ∂ξ + v η ∂v η ∂η ) + 1 ρ ∂p ∂η = g η (3.43)
Now, assume that all the gradients in direction e η are zero, which means that between the particle i and j, the fluid variables are constant in direction e η . Hence, only the effect in direction n ij or e ξ is taken into account when studying the pair of fluid particles i and j, which leads to a quit logical approximation, because the influence in direction e η for the particle i will be compensated when studying the interaction effect between i and the other particles.

With this approximation, (3.43) is simplified and written as

               ∂ρ ∂t + ρ ∂v ξ ∂ξ + v ξ ∂ρ ∂ξ = 0 ∂v ξ ∂t + v ξ ∂v ξ ∂ξ + 1 ρ ∂p ∂ξ = g ξ ∂v η ∂t + v ξ ∂v η ∂ξ = g η (3.44)
As the local speed of sound is defined as

∂p ∂ρ = c 2 (3.45) hence, ∂p ∂ξ = ∂p ∂ρ ∂ρ ∂ξ = c 2 ∂ρ ∂ξ (3.46)
with which (3.44) becomes

               ∂ρ ∂t + ρ ∂v ξ ∂ξ + v ξ ∂ρ ∂ξ = 0 ∂v ξ ∂t + v ξ ∂v ξ ∂ξ + c 2 ρ ∂ρ ∂ξ = g ξ ∂v η ∂t + v ξ ∂v η ∂ξ = g η (3.47)
or in a more concise form

∂W ∂t + A ∂W ∂ξ = Q (3.48)
with,

W =    ρ v ξ v η    , A =    v ξ ρ 0 c 2 /ρ v ξ 0 0 0 v ξ    , and Q =    0 g ξ g η    (3.49)
This quasi-linear system of first order partial differential equations (3.48) is hyperbolic, because its homogeneous part admits wave-like solutions [START_REF] Hirsch | Numerical computation of internal and external flows[END_REF], which means that the matrix A possesses only the real eigenvalues obtained by solving det|A -λI| = 0

       λ 1 = v ξ -c λ 2 = v ξ λ 3 = v ξ + c (3.50)
which correspond to three associated waves in the ξ -t plane, as shown in Figure 3.3, divided into three specific regions: the left region, the right region and a star region denoted by * . W L and W R denote the state of variables for the left and the right regions, respectively. This is a classical Riemann problem where the discontinuity is located at the point O, which has the exact solution. For more details for finding this exact solution, the book [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics[END_REF] of Professor E. F. Toro is highly recommended. Shortly speaking, the wave associated with the eigenvalue λ 2 = v ξ is a contact wave across which we have the continuities

ρ * L = ρ * R = ρ * v * L,ξ = v * R,ξ = v * ξ (3.51)
The other two waves may be shock wave or rarefaction wave, depending on whether it satisfies the entropy condition [START_REF] Chorin | A mathematical introduction to fluid mechanics[END_REF][START_REF] Leveque | Numerical methods for conservation laws[END_REF][START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics[END_REF]. For a shock wave, the Rankine-Hugoniot condition [START_REF] Rankine | On the thermodynamic theory of waves of finite longitudinal disturbances[END_REF] is applied to determine the relation of variables between the left or right region and the star region, whereas for a rarefaction wave, the Generalized Riemann Invariants condition will be used.

As presented in [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics[END_REF], the Generalized Riemann Invariants (GRI) conditions are relations that hold true across the wave structure. The GRI conditions are given as

   dρ + ρ c dv ξ = 0 across λ 1 -wave dρ - ρ c dv ξ = 0 across λ 3 -wave (3.52)
or for the pressure dp + ρcdv ξ = 0 across λ 1 -wave dp -ρcdv ξ = 0 across λ 3 -wave (3.53)

Remark: The speed of wave in Figure 3.3 is not equal to the eigenvalue. Because the system of partial differential equations is non-linear, the eigenvalues are different in the left region and the right region. Hence, one should first detect the type of the wave, then use different methods to determine the speed of the wave. Furthermore, as the investigated flows are at low mach numbers, the speed of flows is markedly smaller than the one of propagation of sound waves. Therefore, in such cases, we may only consider the star region.

Since resolving the Riemann problem analytically for each pair of fluid particles requires an iterating procedure which costs lots of computer sources for the entire fluid domain, and is found not necessary through the practical experience [START_REF] Marongiu | Méthode numérique lagrangienne pour la simulation d'écoulements à surface libre -Application aux turbines Pelton[END_REF]. In this case, the use of some approximate Riemann solvers is reasonable and can save the time of numerical simulation. Here two approximate Riemann solvers will be presented, which are the Primitive Variable Riemann Solver (PVRS) [START_REF] Ivings | On Riemann solvers for compressible liquids[END_REF] and the acoustic Riemann solver, although there exist lots of other approximate Riemann solvers, e.g. the Roe's approximate Riemann solver [START_REF] Roe | Approximate riemann solvers, parameter vectors, and difference schemes[END_REF] and the Osher's approximate Riemann solver [START_REF] Osher | Upwind difference schemes for hyperbolic systems of conservation laws[END_REF].

In the first approximate Riemann solver, PVRS, the matrix A in (3.48) is replaced by a constant matrix

A =    v ξ ρ 0 c 2 /ρ v ξ 0 0 0 v ξ    (3.54)
where

             ρ = 1 2 (ρ L + ρ R ) v ξ = 1 2 (v L,ξ + v R,ξ ) c = 1 2 (c L + c R ) (3.55)
Instead of solving a non-linear Riemann problem, we are going to find the exact solution of a linearized Riemann problem governed by

∂W ∂t + A ∂W ∂ξ = 0 (3.56)
This linear Riemann problem is easy to handle and its solution can be given directly

[78]      ρ * = 1 2 (ρ L + ρ R ) + ρ 2c (v L,ξ -v R,ξ ) v * ξ = 1 2 (v L,ξ + v R,ξ ) + c 2ρ (ρ L -ρ R ) (3.57) where v L,ξ = v i • n ij and v R,ξ = v j • n ij .
The acoustic Riemann solver is based on the linearization of the coefficients in the Generalized Riemann Invariants relations (3.52) or (3.53). In such a procedure, the approximated GRI conditions are

     dρ + ρ L c L dv ξ = 0 across λ 1 -wave dρ - ρ R c R dv ξ = 0 across λ 3 -wave (3.58)
with which one can obtain

     (ρ L -ρ * ) + ρ L c L (v L,ξ -v * ξ ) = 0 (ρ R -ρ * ) - ρ R c R (v R,ξ -v * ξ ) = 0 (3.59)
which gives the solution of the state of fluid variables in the star region

       ρ * = c L + c R + v L,ξ -v R,ξ c L /ρ L + c R /ρ R v * ξ = ρ L -ρ R + ρ L v L,ξ /c L + ρ R v R,ξ /c R ρ L /c L + ρ R /c R (3.60)
then the pressure p * will be calculated by the Tait equation (2.44). Finally, the approximate Riemann solver can offer a sufficiently accurate and rapid solution of the investigated Riemann problem between the two states of fluid Φ i and Φ j . With this solution, the intermediate state can be calculated by

       ρ E,ij = ρ * v E,ij = v * ξ n ij p E,ij = p * (3.61)
which are used to calculate G E,ij in (3.39). This approach proposed by Vila is called the SPH-ALE method, because it combines the Smoothed Particle Hydrodynamics method (SPH) and the Arbitrary Lagrangian-Eulerian (ALE) method. Hence, essentially it is a particle method, but it is not necessary for the particle to be a material fluid volume, which is the case in standard SPH methods. In this approach, the fluid particle can move with an arbitrary velocity v. When v is imposed to be equal to the material velocity v f , the fluid particle is a material volume which is followed as it moves through time and space. This feature allows us to naturally track the interface with other medium, e.g. for free surface flows or interaction with solid boundary. Imposing v = 0, all the fluid particles are fixed in space, which means that an Eulerian mesh is used to discretize the fluid sub-domain. This is similar with the typical Finite Volume (FV) method, which is used to study the flows' properties in some fixed region. The difference is that in this SPH-ALE method the fluid flux between the elemental volumes is calculated with taking into account the interaction effects with all the particles near to the investigated one, whereas in the FV method the flux is obtained by using a certain numerical scheme just at the interface between the elemental volumes. Another important feature of the SPH-ALE method is the use of Riemann solver, which introduces implicitly the artificial viscosity to stabilize and improve the numerical simulation.

Higher order scheme

In the previous subsection, the intermediate status of fluid is obtained by solving approximately a Riemann problem between two fluid particles. This can be classified as first order Godunov-type numerical scheme which was initially introduced by Godunov in 1959 [START_REF] Godunov | A difference scheme for numerical computation of discontinuous solution of hydrodynamic equations[END_REF]. In such type of method, the solution is considered as piece-wise constant over each mesh cell at a fixed time, as shown in Figure 3.4 and the evolution of the flow to the next time step results from the wave interactions originating at the boundaries between adjacent cells [START_REF] Hirsch | Numerical computation of internal and external flows[END_REF].

Through the practical experience of using this first order scheme in space [START_REF] Leduc | Étude physique et numérique de l'écoulement dans un système d'injection de turbine Pelton[END_REF][START_REF] Marongiu | Méthode numérique lagrangienne pour la simulation d'écoulements à surface libre -Application aux turbines Pelton[END_REF], too much artificial viscosity is introduced to the numerical solution, which may interfere the quality of the numerical results. A first possible solution is to use more fluid particles for the spatial discretization. However, this may increase largely the time of calculation. Another solution is the use of higher order numerical schemes to represent the real physical variable fields more accurately.

In 1979 [START_REF] Van Leer | Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov's method[END_REF], van Leer introduced a second order Godunov-type scheme, which is called Monotonic Upstream-centered Scheme for Conservation Laws (MUSCL). In this method, the physical variable is considered to be linearly distributed in each mesh cell, instead of the piece-wise distribution hypothesis for the first order Godunov-type scheme. Hence, the procedure of second order scheme applied in the SPH-ALE method is summarized as follows:

• At the begging of each time-step or sub-time-step, the fluid domain is considered to be discretized into a number of particles;

• When studying the interaction effect between two particles, the mono-dimensional Riemann problem is solved by an approximate Riemann solver between the two status of fluid which are obtained by a linear extrapolation in space;

• The results in the star region of the Riemann problem is then considered as the intermediate status between the two investigated particles. Then G E,ij is calculated by this intermediate status.

Boundary conditions treatment

For solving partial differential equations, the boundary conditions must be given properly, which is often called well-posed, in order to ensure the existence and uniqueness for the analytical solution. Hence the boundary conditions treatment is an essential aspect of the practical application of a numerical scheme.

In the SPH method, the treatment of solid wall boundary condition is a difficulty obstructing the development of this numerical method. As the fluid is considered to be inviscid, the solid wall boundary is treated with the slipping condition which imposes only the equality of velocity in the normal direction of the solid boundary

v f (x k ) -v s (x k ) • n k = 0 (3.62)
where x k denotes the position of the solid wall boundary, and n k is the normal vector.

In the literature, there exist different methods for treating this sort of solid wall boundary condition (3.62), for instance, the repulsive force method [START_REF] Monaghan | Simulating gravity currents with SPH III Boundary forces[END_REF], the fictive particle method [START_REF] Morris | Modeling low Reynolds number incompressible flows using SPH[END_REF], and the ghost particle method [START_REF] Guilcher | Contribution au développement d'une méthode SPH pour la simualtion d'interaction houle-structure[END_REF]. However, these methods appear either inaccurate (the repulsive method), or difficult to handle the complex geometry problems (fictive or ghost particle method). In 2010, Marongiu [START_REF] Marongiu | Free surface flows simulations in Pelton turbines using an hybrid SPH-ALE method[END_REF] introduced an original method for treating the moving solid wall boundary condition. In this method, the boundary effect is obtained by solving a partial Riemann problem at the solid wall, which is coherent with the use of Riemann solver for studying the interaction effect between two fluid particles. Furthermore, this method allows us to easily handle the problems with complex geometries, since we do not need to set extra fictive particles in the solid sub-domain.

Besides, there exist other types of boundary conditions, such as the free surface boundary condition, or the multi-fluids interface boundary condition. While these will not be presented in this work, because it is the solid wall boundary condition that is of interest for the fluid-structure interaction, which is the main subject of this PhD work. More details for these sorts of boundary conditions can be found in [START_REF] Leduc | Étude physique et numérique de l'écoulement dans un système d'injection de turbine Pelton[END_REF][START_REF] Marongiu | Méthode numérique lagrangienne pour la simulation d'écoulements à surface libre -Application aux turbines Pelton[END_REF].

Hyperbolic propagation-dominated system

As presented in Chapter 2, the fluid sub-domain is governed by the Euler equations which possesses the features of hyperbolic equations. This means that the physical problems to be studied are propagation-dominated, but not dominated by diffusion phenomena. The homogeneous part of the governing equations writes

∂U ∂t + A ∂U ∂x + B ∂U ∂y = 0 (3.63)
with,

U =    ρ v x v y    , A =    v x ρ 0 c 2 /ρ v x 0 0 0 v x    and B =    v y ρ 0 0 v y 0 c 2 /ρ 0 v y    (3.64)
This system of equations possesses three real eigenvalues [START_REF] Hirsch | Numerical computation of internal and external flows[END_REF][START_REF] Hirsch | Numerical computation of internal and external flows[END_REF]]

       λ 1 = v • κ -c|κ| λ 2 = v • κ λ 3 = v • κ + c|κ| (3.65)
where κ is the direction vector of the propagation of the wave-like solution, which is usually chosen to be the normal vector n i when studying the boundary condition at Γ i .

At the boundary Γ i , a positive eigenvalue associates with a wave that brings the information from the outside into the inside of the fluid domain, whereas a negative eigenvalue corresponds to a wave bringing the information from the inside to the outside of the fluid domain. Hence, at a certain boundary, if there are totally m = p + n eigenvalues, where p is the number of positive eigenvalues, and n is the one of negative eigenvalues, one should impose p physical boundary conditions, and n numerical boundary conditions.

Imposing a physical boundary condition means imposing the value of certain physical variable or imposing a relationship among several physical variables. For instance, at the entry of the domain Γ 1 in Figure 3.5, the three eigenvalues are

       λ 1 1 = v • n 1 -c λ 1 2 = v • n 1 λ 1 3 = v • n 1 + c (3.66)
where v • n 1 > 0, hence if at the entry the flow is subsonic, there will be two positive eigenvalues λ 2 and λ 3 , and one negative eigenvalue λ 1 . In this case, one can only impose two boundary conditions at this entry, for instance, the fluid density ρ = ρ B1 and the component in direction x of the fluid velocity

v x = v B1
x , or any other combination of two variables among ρ, v x and v y . However, one can not impose all the three physical variables at the entry Γ 1 , since there are only two waves bringing the information from the outside to the inside of the fluid domain. Hence, the third variable that is not imposed should be given by some numerical methods, i.e. the numerical boundary condition. For example, the linear or 2 nd order extrapolation method can give the value of the third variable by using the information coming from the inside of the fluid domain. In contrast, if all the three variables are imposed, this boundary condition may be noncompatible with the governing equations used in the fluid domain.

The method for imposing the numerical boundary conditions is very important and can not be emphasized enough. From the theoretical point of view, as well as through numerical experiments, it can be shown that the choice of the numerical boundary conditions may have a dominating effect on the accuracy, numerical stability and convergence rate of many numerical schemes. For instance, an unconditionally stable scheme may appear to be only conditionally stable in practice, if an incompatible boundary condition is introduced [START_REF] Hirsch | Numerical computation of internal and external flows[END_REF].

Truncation of the support domain of kernel function

Chapter 3.1.1 presents the basis of the SPH method. In (3.13), the term associated with the solid wall boundary is neglected since the kernel function W (x ′x, h) = 0 for x ′ ∈ ∂D, which is true only for the fluid particles sufficiently far from any solid wall boundary. Hence, the semi-discrete equations (3.39) developed previously are not adoptable for the fluid particles which are located near to the solid wall boundary.

A fluid particle is said to be near to the solid wall means that the support domain of its kernel function is truncated by the solid wall, as the fluid particle i shown in Figure 3 For this type of fluid particles, the boundary term in (3.12) can no longer be neglected, because the kernel function W is not equal to zero at the solid wall boundary. Hence the kernel approximation of a gradient function writes

∇ x f (x) ≃< ∇ x f (x) > = D ∇ x ′ f (x ′ )W (x -x ′ , h) dx ′ = ∂D [f (x ′ )W (x -x ′ , h)]n dx ′ - D f (x ′ )∇ x ′ W (x -x ′ , h) dx ′ (3.67)
Using the particle approximation gives the gradient for the fluid particle i

∇ x f (x)| x=x i ≃ < ∇ x f (x) >| x=x i ≃ k∈∂D i f (x k )W (x i -x k , h)n k s k - j∈D i f (x j )∇ j W (x i -x j , h)ω j = k∈∂D i f k W ik n k s k - j∈D i f j ∇ j W ij ω j = k∈∂D i f k W ik n k s k + j∈D i f j ∇ i W ij ω j (3.68)
where k denotes the discrete surface element, n k is the unit normal vector pointing from the fluid to the solid, as shown in Figure 3.6, and s k is the surface area of this element.

Remark: In this work, the "surface" element is assumed to be either a point (1-D) or a straight line (2-D), or a plane surface (3-D). Hence the normal vector

n k is constant inside a "surface" element k. When f (x) is the flux tensor in (3.29), [∇ x • F] i rewrites [∇ x • F] i ≃ k∈∂D i W ik F k • n k s k + j∈D i F j • ∇ i W ij ω j (3.69)
Now, applying the SPH approximation (3.68) to calculate the inserted term (3.36) in (3.35), which is used to keep the conservation property

Ω F i • ∇ x H dΩ = F i • Ω ∇ x H dΩ ≃ F i • ω i k∈∂D i H k W ik n k s k + ω i j∈D i H j ∇ i W ij ω j = ω i k∈∂D i W ik F i • n k s k + ω i j∈D i F i • ∇ i W ij ω j
d(ω i Φ i ) dt + ω i k∈∂D i W ik (F k + F i ) • n k s k + ω i j∈D i (F j + F i ) • ∇ i W ij ω j = S i ω i (3.
71) which has an additional term associated with the effect of solid boundary, comparing with (3.37). This equation has the same form with one used in [START_REF] De Leffe | Modélisation d'écoulements visqueux par méthode SPH en vue d'application à l'hydrodynamique navale[END_REF], although the derivation procedure is different.

As presented in Chapter 3.1.4, Vila [START_REF] Vila | On particle weighted methods and smooth particle hydrodynamics[END_REF] proposed to approximate (F j + F i ) by 2G E,ij (Φ i , Φ j ) which is calculated with an intermediate status of fluid variable obtained by solving a Riemann problem between the two fluid particles i and j. In 2010 [START_REF] Marongiu | Free surface flows simulations in Pelton turbines using an hybrid SPH-ALE method[END_REF], Marongiu introduced a method for treating the moving solid wall boundary condition. In this method, (F k + F i ) is replaced by 2G E,ik (Φ i , Φ k ) which is calculated with an intermediate status obtained by solving a partial Riemann problem. Hence, (3.71) is rewritten as

d(ω i Φ i ) dt + ω i k∈∂D i W ik 2G E,ik • n k s k + ω i j∈D i 2G E,ij • ∇ i W ij ω j = S i ω i (3.72)
Taking Φ i = ρ i and Φ i = ρ i v i , one can obtain the semi-discrete equations of conservation laws for the mass and the linear momentum

                             d(ρ i ω i ) dt + ω i k∈∂D i W ik 2ρ E,ik (v E,ik -v ik ) • n k s k + ω i j∈D i 2ρ E,ij (v E,ij -v ij ) • ∇ i W ij ω j = 0 d(ρ i ω i v i ) dt + ω i k∈∂D i W ik 2 [ρ E,ik v E,ik ⊗ (v E,ik -v ik ) + p E,ik I] • n k s k + ω i j∈D i 2 [ρ E,ij v E,ij ⊗ (v E,ij -v ij ) + p E,ij I] • ∇ i W ij ω j = ρ i g i ω i (3.73)
where ρ E,ik , v E,ik and p E,ik denote the intermediate fluid status obtained by solving a partial Riemann problem at the solid wall boundary. In this work, the Lagrangian mode is chosen to be used for coupling with the struc-ture. In such mode

v i = v i and v ij = 1 2 (v i + v j ). Moreover, we impose v E,ik = v ik which leads to                    d(ρ i ω i ) dt + ω i j∈D i 2ρ E,ij (v E,ij -v ij ) • ∇ i W ij ω j = 0 d(ρ i ω i v i ) dt + ω i k∈∂D i W ik 2p E,ik n k s k + ω i j∈D i 2 [ρ E,ij v E,ij ⊗ (v E,ij -v ij ) + p E,ij I] • ∇ i W ij ω j = ρ i g i ω i (3.74)
One of the advantages of imposing v E,ik = v ik is that with the equation of conservation of mass in (3.74) we can enforce no mass transfer at the solid wall boundary, which is a very important physical feature for the numerical method. Finally, the only value to be solved at the solid wall boundary is the intermediate status of pressure p E,ik .

Remark: In Vila's formula [START_REF] Vila | On particle weighted methods and smooth particle hydrodynamics[END_REF], the condition v E,ij = v ij is not necessary, which means that even the Lagrangian mode is used, the mass transfer may occur between two fluid particles. As Vila commented, although we loose this nice property of the method, we keep global conservation of mass and momentum.

Partial Riemann problem

As mentioned in the previous part, the intermediate fluid status ρ E,ik , v E,ik and p E,ik are obtained by solving a partial Riemann problem. This approach was introduced by Marongiu [START_REF] Marongiu | Méthode numérique lagrangienne pour la simulation d'écoulements à surface libre -Application aux turbines Pelton[END_REF][START_REF] Marongiu | Free surface flows simulations in Pelton turbines using an hybrid SPH-ALE method[END_REF] in the SPH method for treating the moving solid boundary condition.

Consider a solid surface element k truncating the support domain the fluid particle i, as shown in Figure 3.7, n k is the normal vector of this surface element. Because the fluid is considered to be inviscid, the slipping condition (3.62) is imposed. In other words, when treating the solid boundary condition, only the effect in the normal direction n k should be taken into account. For this reason, the partial Riemann problem will be resolved in the ξ -t plane illustrated in Figure 3.7.

For the classical Riemann problem, such as in Section 3.1.5, the ξ-t plane is divided into three regions by two waves associated with the two eigenvalues v ξ -c and v ξ + c. Then the Generalized Riemann Invariants conditions or the Rankine-Hugoniot relation will be used to relate the status in the star region with the one in the left or right region in order to determine the status in this star region, which will be used as an intermediate status. However, in the case presented in Figure 3.7, we can have only one relation across the wave associated with the eigenvalue v ξ -c, because there is no fluid particle inside of the solid wall. Hence we need another relation between the variables ρ and v ξ . This is called the partial Riemann problem [START_REF] Dubois | Partial riemann problem, boundary conditions, and gas dynamics. Absorbing Boundaries and Layers, Domain Decomposition Methods: Applications to Large Scale Computations[END_REF], in which one needs to impose a relation or manifold M(ρ * , v * ξ ) = 0 to determine the status in the star region. In fact, this manifold plays a role of the physical boundary condition for the hyperbolic system of equations. As presented in [START_REF] Marongiu | Free surface flows simulations in Pelton turbines using an hybrid SPH-ALE method[END_REF], across the wave associated with v ξ -c, a relation between ρ and v ξ can be found and expressed as a curve in the ρ-v ξ plane, as illustrated in Figure 3

p E,ik = p * = p i -ρ i c i (v k -v i ) • n k (3.79)
with which we can calculate the boundary term

ω i k∈∂D i W ik 2p E,ik n k s k in (3.74) in order
to take into account the effect of moving or fixed solid wall boundary.

Fluid pressure at solid wall

In 2007, Marongiu [START_REF] Marongiu | Méthode numérique lagrangienne pour la simulation d'écoulements à surface libre -Application aux turbines Pelton[END_REF] introduced a method for calculating the fluid pressure at the solid wall boundary to estimate the pressure field on the impacted surface of the blade of Pelton turbine. When Φ i = ρ i v i , (3.72) represents the conservation equation of the linear momentum for the fluid particle i which is located near the solid wall boundary. It is obvious that the linear momentum of the fluid particle i varies in time due to three types of sources • The solid wall boundary force

-ω i k∈∂D i W ik 2G E,ik • n k s k (3.80)
• The interaction effect between the fluid particles

-ω i j∈D i 2G E,ij • ∇ i W ij ω j (3.81)
• The body force (gravity) effect ρ i g i ω i

As shown previously, in Lagrangian mode G E,ik = p E,ik I by imposing v E,ik = v ij . Hence, for the fluid particle i, the force exerted by the solid surface element k is

F k→i = -ω i W ik 2p E,ik n k s k (3.82)
Inversely, the Newton's third law tells us that the force applied to the solid surface element k from the fluid particle i is

F i→k = -F k→i = ω i W ik 2p E,ik n k s k (3.83)
then the force exerted to solid surface element k from the whole fluid domain Ω f can be calculated by

F Ω f →k = i∈D k F i→k = i∈D k ω i W ik 2p E,ik n k s k (3.84)
where i ∈ D k denotes all the fluid particles which are sufficiently near to the solid surface element k, i.e. the support domain of the kernel function is truncated by the solid wall boundary.

Given the definition of the fluid pressure p k at the solid surface element k

F Ω f →k = p k n k s k (3.85)
we can obtain

p k n k s k = i∈D k ω i W ik 2p E,ik n k s k (3.86)
then the fluid pressure p k can be calculated by

p k = i∈D k ω i W ik 2p E,ik (3.87)

System of semi-discrete equations for fluid

In Section 3.2.2, we obtained the semi-discrete equations (3.72) for the conservation of mass and linear momentum. The temporal evolution equation of the position of fluid particles is given as

dx i dt = v i (3.88)
where v i denotes the arbitrary velocity of the control volume Ω i . In Lagrangian mode, v i = v i , whereas in Eulerian mode, v i = 0.

In addition, another equation is needed to describe the evolution of the volume of Ω i . For doing this, let's take f (x, t) = 1, ∀x ∈ R d , ∀t ∈ R + , in the Reynolds Transport Theorem (2.25), we obtain

∂ω i ∂t χ = Γ i v(x, t) • n dx = Ω i ∇ x • v dx = ω i [∇ x • v] i (3.89)
where ω i denotes the volume of the fluid particle Ω i , and

[∇ x • v] i is the mean value of the divergence of v over Ω i .
Using the kernel and particle approximations (3.68) gives the semi-discrete equation

for ω i dω i dt = ω i k∈∂D i W ik v k • n k s k + ω i j∈D i v k • ∇ i W ij ω j (3.90)
To ensure a zero order consistency for the approximation of ∇ x • v, we insert an additional term into (3.89)

∂ω i ∂t χ = Ω i (∇ x • v -v i • ∇ x H) dΩ = ω i [∇ x • v -v i • ∇ x H] i (3.91)
where

H(x) = 1 ∀x ∈ R d . Since ∇ x H = 0, (3.91) is mathematically identical with (3.89).
Using the same procedure (3.70) to discretize the inserted term, we can obtain the semi-discrete form of (3.91)

dω i dt = ω i k∈∂D i W ik ( v k -v i ) • n k s k + ω i j∈D i ∇ i W ij • ( v j -v i )ω j (3.92)
which describes the evolution in time of the volume ω i of the fluid particle. Gathering (3.73), (3.88) and (3.92), we get the whole system of semi-discrete equations for the fluid sub-domain

                                               dx i dt = v i dω i dt = ω i k∈∂D i W ik ( v k -v i ) • n k s k + ω i j∈D i ∇ i W ij • ( v j -v i )ω j d(ρ i ω i ) dt + ω i k∈∂D i W ik 2ρ E,ik (v E,ik -v ik ) • n k s k + ω i j∈D i 2ρ E,ij (v E,ij -v ij ) • ∇ i W ij ω j = 0 d(ρ i ω i v i ) dt + ω i k∈∂D i W ik 2 [ρ E,ik v E,ik ⊗ (v E,ik -v ik ) + p E,ik I] • n k s k + ω i j∈D i 2 [ρ E,ij v E,ij ⊗ (v E,ij -v ij ) + p E,ij I] • ∇ i W ij ω j = ρ i g i ω i (3.93)
Choosing the Lagrangian mode, (3.93) becomes

                                       dx i dt = v i dω i dt = ω i k∈∂D i W ik (v k -v i ) • n k s k + ω i j∈D i ∇ i W ij • (v j -v i )ω j d(ρ i ω i ) dt + ω i j∈D i 2ρ E,ij (v E,ij -v ij ) • ∇ i W ij ω j = 0 d(ρ i ω i v i ) dt + ω i k∈∂D i W ik 2p E,ik n k s k + ω i j∈D i 2 [ρ E,ij v E,ij ⊗ (v E,ij -v ij ) + p E,ij I] • ∇ i W ij ω j = ρ i g i ω i
(3.94)

Explicit Runge-Kutta time integrator

The SPH method is used only for the spatial discretization of the fluid sub-domain which is governed by the partial differential equations. The system of equations (3.94) is called semi-discrete because the time domain is also needed to be discretized, or these temporal differential equations should be integrated in time by using certain numerical time integrators. In this work, the 2 nd order explicit Runge-Kutta scheme is used as the time integrator. After the spatial discretization, the equations in (3.94) are all ordinary temporal differential equations, which can also be written as

dU dt = H(U ) (3.95)
where U = U (t) depending only on time t. In the system (3.94) for the fluid particle i, U may be the position x i , the volume ω i , the mass ρ i ω i or the momentum

ρ i ω i v i . H(U )
is a function of U , which represents the spatial discretization part. Assume that at the instant t = t n , all the fluid variables U n are already known, the objective is to calculate the state of variables U n+1 at the next instant t = t n+1 , as shown in Figure 3.9.

The 2 nd order explicit Runge-Kutta scheme is often called the two-stage Runge-Kutta scheme, because there are two stages for going to t n+1 from t n . Figure 3.9: Advancing in time from the instant t n to t n+1 by using time integrator.

• Stage 1: Calculate the state of variables U n+1/2 at the instant t n+1/2

U n+1/2 = U n + ∆t 2 H(U n ) (3.96)
where U n+1/2 denotes the state of variables at the instant t n+1/2 , which is obtained by using a first order time integrator.

• Stage 2: Calculate the state of variables U n+1 at the instant t n+1

U n+1 = U n + ∆tH(U n+1/2 ) (3.97)
where the function H(U n+1/2 ) is calculated with the variables U n+1/2 obtained in the first stage.

Because one should calculate first the state of variables at the instant t n+1/2 , this method is also called the mid-point version of the 2 nd order Runge-Kutta scheme. However, the 2 nd order Runge-Kutta scheme does not have a unique form, for example, the Heun's method is also classified as a 2 nd order Runge-Kutta scheme. In this work, it is the mid-point version that is used for time integration of the fluid equations and will be coupled with the time integrator used in the solid sub-domain.

As an explicit time integrator is used, the Courant-Friedriech-Levy (CFL) condition needs to be satisfied in order to have a stable numerical simulation. In this work, the applied CFL condition is the one given in [START_REF] Leduc | Étude physique et numérique de l'écoulement dans un système d'injection de turbine Pelton[END_REF] ∆t

K CFL min i h i c i + 2|v i | (3.98)
where the coefficient K CFL ∈ (0, 1) is in function of the used time integrator, which is usually determined through the practical numerical experience.

Chapter 4

Finite Element numerical method for solid

The mechanical phenomena in the solid sub-domain are governed by some partial differential equations (2.45) or (2.47). In general, solving analytically these equations is difficult or even impossible. Hence, numerical methods are often used to obtain approximate solutions of these governing equations in engineering as well as in academic research. This chapter starts with a brief introduction of a widely used numerical method, the Finite Element (FE) method [START_REF] Belytschko | Nonlinear Finite Elements for Continua and Structures[END_REF][START_REF] Fish | A first course in finite elements[END_REF]36], in solid mechanics. Secondly, the treatments of different types of boundary conditions are presented. After the spatial discretization by the FE method, the third part gives the semi-discrete system of equations in matrix form. Finally, the used time integrator is presented at the end of this chapter.

Finite Element method

The Finite Element (FE) method is a numerical approach by which the partial differential governing equations can be solved approximately [START_REF] Fish | A first course in finite elements[END_REF]. Comparing to other numerical methods, e.g. the Finite Difference (FD) method, one of the advantages of using the FE method is the flexibility for complicated geometries.

The basic idea of the FE method in solid mechanics is to divide the body Ω s into finite elements

Ω s = Ns e Ω e (4.1)
where N s denotes the number of the elements, Ω e is one of the elemental volumes. All the elements are connected by the nodes, and this is called the computational mesh or mesh. The FE method can give the approximate solutions of the governing equations at these mesh points.

For linear problems, the solution is determined by solving a system of linear equations for which the number of unknowns is equal to the number of nodes. For the FE method in solid mechanics, there are three major sources of non-linearity:

• Large deformation and displacement (non-linear geometry)

• Non-linear constitutive model (non-linear material)

• Contact problems

In this work only the non-linear geometry and non-linear material will be considered. The solid-solid contact problems are not investigated in this PhD work.

Weak form for updated Lagrangian formulation

As presented in Chapter 2, the FE discretizations are commonly classified as Updated Lagrangian Formulations (ULF) and Total Lagrangian Formulations (TLF) [START_REF] Belytschko | Nonlinear Finite Elements for Continua and Structures[END_REF]. Notice that both the two formulations are with Lagrangian meshes, which means that the nodes and elements move with the material, and the boundaries and interfaces remain coincident with element edges. Furthermore, the constitutive equations are always evaluated at the same material points, which is advantageous for history-dependent materials. More concretely, the time derivatives in the governing equations (2.45) or (2.47) are both carried out holding the material coordinates constant. The differences between the two formulations are that in ULF the spatial derivative is with respect to the spatial coordinates and the corresponding weak form is obtained by integration over the current configuration of the solid body, whereas in TLF the spatial derivative is with respect to the material coordinates and its weak form is given by integration over the initial configuration. Because the ULF and TLF are mathematically equivalent, only the weak form for the ULF is presented in this section.

• Velocity boundary condition v(X, t) = v b (X, t) on Γ v (4.2)
• Traction boundary condition

σ(X, t) • n b (X, t) = t b (X, t) on Γ t (4.3)
where v(X, t) denotes the velocity vector expressed with the material coordinates X, v b is the imposed velocity at the boundary Γ v . σ(X, t) represents the stress tensor, n b the normal vector of the boundary Γ t , and t b (X, t) the traction or force exerted on Γ t . The relation between Γ v and Γ t is described as

Γ v ∪ Γ t = Γ s Γ v ∩ Γ t = ∅ (4.4)
where Γ s denotes the boundary of the whole solid domain Ω s . The weak form is obtained by multiplying the momentum equation by a test function and then integrating over the current domain, which gives

Ωs δv • ρ dv dt -∇ x • σ -ρb dΩ = 0 (4.5)
where δv(X) is the test function satisfying

δv(X) ∈ U 0 U 0 = δv(X) | δv ∈ C 0 (X), δv = 0 on Γ v (4.6)
and the condition on the trial function v(X, t) is

v(X, t) ∈ U U = v(X, t) | v ∈ C 0 (X), v = v b on Γ v (4.7)
where a C 0 function is continuous in space, but its first order derivative is only piecewise differentiable.

Applying the traction boundary condition (4.3) and Gauss's theorem, (4.5) can be written as [START_REF] Belytschko | Nonlinear Finite Elements for Continua and Structures[END_REF] Ωs

(∇ x δv) ⊤ : σ dΩ - Ωs δv • ρb dΩ - Γt δv • t b dΓ + Ωs δv • ρ v dΩ = 0 (4.8)
where (∇ x δv) ⊤ is the spatial gradient tensor of the virtual velocity, which is denoted by δL and can be calculated with the definition of the left gradient operator (2.16) where δD is called the virtual strain rate, which is the symmetric part of δL, and the spin δW is the skew-symmetric part of δL. Because the Cauchy stress tensor σ is symmetric, i.e. σ ⊤ = σ, then we have (∇ x δv) ⊤ : σ = δL : σ = (δD + δW) : σ = δD : σ (4.12)

δL = (∇ x δv) ⊤ =          ∂δv x ∂x
since δW : σ = 0.

The procedure for obtaining the weak form is also called the virtual power principle, in which the test function δv is called the virtual velocity. In this principle, (4.5) or (4.8) implies that for an arbitrary virtual velocity δv, the virtual power δP tot of the whole system is equal to zero. In fact, in (4.8) each part possesses a physical sense: This weak form for Updated Lagrangian Formulation will be discretized in next section.

Finite Element discretization

In the FE method, the motion x(X, t) is approximated by

x(X, t) = n I N I (X)x I (t) (4.17)
where I denotes the identification number of node, and n is the total numbers of the nodes. N I (X) is a scalar function of X, which is called the shape function and x I is the position vector of the node I. With the Einstein's notation method, (4.17) can be simply written as x(X, t) = N I (X)x I (t) (4.18)

Using the same approximation methodology to the displacement field, one obtains

u(X, t) = N I (X)u I (t) (4.19)
where u I (t) = x I (t) -X I with X I being the initial position of the node I.

The velocity is obtained by taking the material time derivative of the displacement v(X, t) = ∂u(X, t) ∂t

X = ∂ [N I (X)u I (t)] ∂t X = N I (X) uI (t) (4.20)
Similarly, the acceleration field is approximated by a(X, t) = ∂v(X, t) ∂t

X = v(X, t) = N I (X)ü I (t) (4.21)
The test function or virtual velocity δv(X) is not a function of time, hence we approximate the test function as

δv(X) = N I (X)δv I (4.22)
with which δL can be written as

δL = ∂δv ∂x = ∂[N I (X)δv I ] ∂x =          ∂N I ∂x
         = δv I ⊗ ∇ x N I (X) (4.23) or in matrix form δL = δv I [∇ x N I (X)] ⊤ (4.24)
with,

δv =    δv x I δv y I δv z I    and ∇ x N I (X) =         ∂N I ∂x ∂N I ∂y ∂N I ∂z         = B I (4.25)
where δv x I , δv y I and δv z I are the components of the virtual velocity δv I of the node I in the three directions.

Remark: In this thesis, a variable or an equation may be expressed either in tensor

form or in matrix form, which are essentially equivalent. When the operation is carried out with an operator like "•" or "⊗" etc., it means that we are using the tensor form. Whereas, if in the equation it has the transpose symbol "⊤", it is in the matrix form.

With the approximation methods presented previously, the three types of virtual power can be calculated as

                 δP int = δv ⊤ I Ωs σB I dΩ = δv ⊤ I f int I δP ext = δv ⊤ I Ωs N I ρb dΩ + Γt N I t b dΓ = δv ⊤ I f ext I δP kin = δv ⊤ I Ωs N I ρN J dΩ vJ = δv ⊤ I M IJ vJ
M n1 I M n2 I • • • M nI I • • • M nn I             , f int =             f int 1 f int 2 . . . f int I . . . f int n             , f ext =             f ext 1 f ext 2 . . . f ext I . . . f ext n            
(4.32) where I denotes the unit matrix of the dimension of the problems, and v is the velocity matrix of dimension [nd × 1] with d being the number of dimensions of the investigated problems and n the total number of nodes

v s = v 1 v 2 . . . v I . . . v n ⊤ (4.33)
and u is the nodal displacement matrix

u s = u 1 u 2 . . . u I . . . u n ⊤ (4.34)

Linear elastic constitutive equations

The semi-discrete equations (4.31) is applicable to linear or non-linear cases, since the methods for calculating the internal and external nodal forces are not specified. For instance, when calculating f int , the material property may be linear or non-linear, and the deformation of the solid may be little or large.

A constitutive model gives the relation between the stress and the strain, which is used to close the system of equations for the solid sub-domain. For measuring the stress and strain, various choices may be applied, e.g. the Green strain E, the strain rate D, the Cauchy stress σ, and the second Piola-Kirchhoff tress S etc.

In this part, we only give the example of linear elastic material with small strain. The constitutive model is used to relate the Cauchy stress tensor σ with the Almansi strain tensor ε

σ = C : ε (4.35)
where C is a fourth order tensor describing the linear elasticity, and C is independent with ε. The Cauchy's strain tensor is calculated by

ε = 1 2 ∇ x u + (∇ x u) ⊤ (4.36)
The Voigt notation method is often used for constitutive equations. In this method, (4. In this linear elastic case, the 3-D Hook's law for isotropic material gives where u I = [u x I u y I u z I ] ⊤ denoting the nodal displacement for the node I. With (4.37), (4.43) and (4.45), the virtual internal power can be calculated by

[C] = E (1 + ν)(1 -2ν)              1 -ν ν ν 0 0 0 ν 1 -ν ν 0 0 0 ν ν 1 -ν 0 0 0 0 0 0 1 -2ν 2 0 0 0 0 0 0 1 -2ν 2 0 0 0 0 0 0 1 -2ν 2              (4.
B I =                  ∂N I ∂x 0 0 0 ∂N I ∂y 0 0 0 ∂N I ∂z ∂N I ∂y ∂N I ∂x 0 0 ∂N I ∂z ∂N I ∂y ∂N I ∂z 0 ∂N I ∂x                
δP int = Ωs (B I δv I ) ⊤ [C] B J u J dΩ = δv ⊤ I Ωs B ⊤ I [C] B J dΩ u J = δv ⊤ I f int I (4.46)
Hence the internal nodal force f int I for the node I is

f int I = Ωs B ⊤ I [C] B J dΩ u J = K IJ u J (4.47)
with which the momentum equation for the node I (4.29) becomes

M IJ üJ + K IJ u J = f ext I (4.48) 65
then the semi-discrete equations for all the nodes in matrix form are

Mü s + Ku s = f ext (4.49)
where Ku s = f int , and K is called the stiffness matrix, with

K IJ = Ωs B ⊤ I [C] B J dΩ (4.50)
Notice that in linear elastic cases with small deformation, the stiffness matrix K can be considered to be constant. The mass matrix is always the same, even for the non-linear cases with large deformation, because the mass matrix is calculated by

M IJ = Ωs N I (X)ρ(x, t)N J (X) dx = Ω 0 s N I (X)ρ(X, t)N J (X)J(X, t) dX = Ω 0 s N I (X)ρ 0 (X)N J (X) dX (4.51)
where ρ 0 (X) denotes the initial density at the material point X, and J(X, t) is the determinant of the gradient matrix F = ∂x ∂X ⊤ . As presented in (2.47), the mass conservation gives ρ 0 (X) = ρ(X, t)J(X, t).

Remark:

The shape function N I (X) is always expressed with the material coordinates X. Sometimes we see that it appears in the integral over the current configuration Ω(x).

In such cases, the shape function is expressed by N I Φ -1 (x, t) with Φ -1 (x, t) is the mapping from x to X. However, this inverse mapping is never carried out explicitly in the practical calculations.

Boundary conditions treatment

In the Finite Element method, the boundary conditions are usually of two types: the essential boundary condition and the nature boundary condition corresponding to the velocity boundary condition and the force boundary condition, respectively, as shown in Figure 4.1.

In the essential boundary conditions, the value of the variable is to be imposed in order to solve the partial differential equations completely. This sort of condition is also called the Dirichlet-type boundary condition.

As for the nature boundary conditions, they are often associated with the higher order derivatives of the variables. This kind of boundary condition is not necessary for solving the differential equations completely.

Essential boundary conditions

In the essential boundary conditions or the so called velocity boundary conditions, we impose the values of nodal velocity at the boundary Γ v . For instance, I = b is a node on the boundary Γ v , and the value of the imposed nodal velocity is

v b = [v x b v y b v z b ] ⊤ .
Usually, a linear system of equations is to be resolved to calculate the velocity field

Av = B (4.52)
where A is the coefficient matrix of dimension [nd×nd] and B is the source term matrix which is of dimension [nd × 1]. A and B are both considered to be known. The aim is to resolve the unknown matrix v of dimension [nd × 1].

In this work, the essential boundary conditions are treated by the elimination method. In this method, if the velocity of the node I = b is imposed, the matrix A and B should be modified:

• Modification to the matrix A -Keep the value of A ii with i = d(b -1) + 1, • • • , d(b -1) + d;
-Set all the other elements in the line i or the column i to be zero.

• Modification to the matrix

B -Set B i = A ii v i with v i = v x b , v y b , or v z b .
Other methods can also be applied to treat the Dirichlet-type boundary condition, e.g. the Lagrange Multiplier method, in which the velocity boundary condition is considered to be a constraint condition when minimizing the power or the Lagrangian of the whole system. Comparing to the elimination method, the Lagrange Multiplier method requires to resolve additional equations for the Lagrange multipliers λ. In such approach, the system of equations to be resolved is often written as

Av + C ⊤ Λ = B Cv = D (4.53)
where C is a matrix of dimension [nm × nd] with m being the number of nodes at the velocity boundary Γ v . Hence, Cv = D represents the essential boundary conditions. The Lagrange multipliers matrix Λ contains the additional unknowns to be solved which are related to the force or pressure applied to the solid domain at this boundary Γ v . The resolution procedure of this system of equation (4.53) is as follows:

• Calculate the condensed matrix H H = C[A] -1 C ⊤ (4.54)
• Calculate the Lagrange multipliers matrix

Λ Λ = [H] -1 C[A] -1 B -D (4.55)
• Calculate the velocity field

v = [A] -1 (B -C ⊤ Λ) (4.56)

Natural boundary conditions

The nature boundary conditions are also called the force boundary conditions. We impose the external forces on the boundary Γ t as shown in Figure 4.1.

In this work, the external forces are classified into two categories: (i) f ext,s denotes the external nodal forces applied directly on the solid, which is often known at each instant of time; (ii) f ext,f represents the external nodal forces associated with the effect of the fluid-structure interaction, for which the exerted force can not be explicitly given.

Hence the external nodal forces in (4.31) is calculated as

f ext = f ext,s + f ext,f (4.57)

Newmark time integrator

Notice that the velocity matrix v in (4.33) and the displacement matrix u in (4.34) contain only the time-depending values, since the governing equations are already discretized by the FE method. (4.31) is the semi-discrete system of equations for the solid sub-domain, which will be resolved in time with the Newmark time integrator [START_REF] Newmark | A method of computation for structural dynamics[END_REF].

Given the linear dynamic equilibrium equation ( 4.49) at the instant

t = t n+1 Mü s (t n+1 ) + Ku s (t n+1 ) = f ext (t n+1 ) (4.58) or Ma n+1 s + Ku n+1 s = f n+1 ext (4.59)
where all the variables at the instant t = t n+1 are considered as unknowns and to be resolved with the time integrator. The general form of the Newmark scheme is

   u n+1 s = u n s + ∆tv n s + ∆t 2 2 (1 -2β)a n s + 2βa n+1 s v n+1 s = v n s + ∆t (1 -γ)a n s + γa n+1 s (4.60)
where u n s denotes the displacement field at the instant t = t n , v n s the velocity field and a n s the acceleration field, which are all already known. ∆t = t n+1 -t n representing the time step. β and γ are the coefficients of the Newmark scheme. Different combinations of β and γ give the different properties of the numerical scheme [START_REF] Bathe | Numerical Methods in Finite Element Analysis[END_REF]:

• Unconditionally unstable, if γ < 1 2 (4.61) • Unconditionally stable, if 1 2 γ 2β (4.62) • Conditionally stable, if γ 1 2 , β < γ 2 , ∆t Ω crit ω max (4.63)
where Ω crit for undamped system is

Ω crit = γ 2 -β -1/2 (4.64)
and ω max corresponds to the highest natural frequency of

K -ω 2 M = 0 (4.65)
For example, taking β = 0.25 and γ = 0.5 is equivalent to assuming that the acceleration has a mean value during one time step, hence this is called the mean acceleration Newmark scheme. An explicit central difference scheme is obtained if we set β = 0 and γ = 0.5. These two schemes are both second order accurate in time O(∆t 2 ). If γ > 0.5, the algorithmic damping effect will be introduced, and the numerical scheme becomes first order accurate in time O(∆t).

Chapter 5

One coupling strategy for FSI Previously, Chapter 3 and Chapter 4 present the numerical methods used for the fluid and solid sub-domains, respectively. The SPH-ALE method is used for the fluid subdomain and the Finite Element method is applied for the solid sub-domain. The contribution of this PhD work consists in proposing an energy conserving approach for coupling two different numerical methods for the simulation of fluid-structure interaction, which can preserve the numerical stability as well as the minimal order of accuracy of the used time integrators in all the sub-domains.

This chapter presents the main contribution of this PhD work, beginning with an introduction of the physical conditions at the fluid structure interface. Then an energy conserving approach is introduced for coupling of sub-domains simulated by means of different discretization methods and different time integrators. As a monolithic procedure, this proposed method requires to resolve a coupled system of equations for updating the status of variables at the interface, which is shown next. Finally, the coupling algorithm is presented to give a survey of the proposed coupling method.

Interface conditions

In continuum mechanics, at the interface of two different mediums, the continuity of physical variables should be imposed, according to different physical problems. In this work, only the kinematic continuity and dynamic equilibrium are considered, which leads to the continuity conditions on the velocity and force at the interface.

For a viscous fluid flow interacting with a deformable structure, the continuity conditions at the interface are usually given as

∀x ∈ Γ I v f (x, t) = v s (x, t) σ f (x, t) • n f (x, t) = -σ s (x, t) • n s (x, t) (5.1)
where Γ I denotes the fluid structure interface, n f (x, t) and n s (x, t) are the normal vectors pointing-out of each sub-domain at the interface, with n f (x, t) = -n s (x, t) as shown in Figure 5.1.

For an inviscid fluid, the continuities of velocity and force are not enforced in the tangential direction, hence the slipping condition or non-penetrating condition is used at the interface

∀x ∈ Γ I n f (x, t) • v f (x, t) = -n s (x, t) • v s (x, t) n f (x, t) • [σ f (x, t) • n f (x, t)] = n s (x, t) • [σ s (x, t) • n s (x, t)] (5.2)
which imposes only the continuities in the normal direction.

Because in this work the fluid is considered to be inviscid, we will apply the interface condition (5.2) for the numerical simulation of fluid-structure interaction. As presented in Chapter 2 by (2.42), the stress tensor in the fluid sub-domain can be rewritten as σ f = -p f I. Because of the slipping condition, we can similarly write the solid stress tensor at the interface as σ s = -p s I with p s being the pressure applied on the solid. Hence the condition at the interface Γ I becomes

n f • v f = -n s • v s p f = p s (5.3)
where -p s = n s • t b with t b being the surface force in traction boundary condition (4.3) for the solid structure.

Energy conservation method

In 1997, Farhat et al. [23] show that if we use a hybrid formulation, such as the Schur's dual formulation, to connect the sub-domains and if the link energy is zero, the nu- merical stability and the accuracy of the Newmark scheme in each sub-domain are not affected by the interfaces.

In 2002, Combescure et al. [START_REF] Combescure | A numerical scheme to couple subdomains with different time-steps for predominantly linear transient analysis[END_REF] proposed to use an energy formulation [START_REF] Hughes | Nonlinear Finite Element Analysis[END_REF] to prove that the numerical stability of the global problem depends on the continuity conditions prescribed at the interface. They show that a continuity constraint on velocities at the interface is the best choice to couple sub-domains with different time-steps and different Newmark time integrators.

Based on the work of Combescure, Mahjoubi et al. [START_REF] Mahjoubi | Coupling subdomains with heterogeneous time integrators and incompatible time steps[END_REF]49] proposed a monolithic energy conserving method to couple heterogeneous time integrators with incompatible time steps in structural dynamics. The key point of this proposed method is to ensure the interface conditions, velocity continuity and dynamic equilibrium, in a weak sense in time. This means that instead of verifying the interface conditions at each instant in time, the interface or link energy is ensured to be zero during a time step, and thusly for the whole period of the numerical simulation.

Time lag for staggered methods

As presented in the chapter of introduction, in a monolithic procedure, the interaction of the fluid and the structure at the mutual interface is treated synchronously, which means that at the interface the states of fluid and structure are obtained at the same time. This leads to the need of resolving the equations of different sub-domains in one simulation software, or modifying sweepingly the two existing calculation codes. However, the advantage of monolithic approach is the preservation of numerical stability and accuracy.

Partitioned procedure can overcome this limit, because the fluid and structure equations are alternatively integrated in time and the interface conditions are enforced asyn-chronously [START_REF] Michler | A monolithic approach to fluid-structure interaction[END_REF]. A general survey of partitioned procedures is given in [26] by Felippa et al. In 1995, Piperno et al. [START_REF] Piperno | Partitioned procedures for the transient solution of coupled aeroelastic problems -Part I: Model problem, theory and two-dimensional application[END_REF] presented several partitioned schemes for fluidstructure interaction with Arbitrary-Lagrangian-Eulerian (ALE) formulation.

Assume that the time domain is discretized into time steps, and at the instant t n all the physical variables of the fluid and the structure are already known, the objective is to calculate the status of variables at the next instant t n+1 . In such case, the most typical and basic partitioned procedure is carried out on the following sequential process [69]: The procedure of basic partitioned coupling method. F n and S n denote the status of the fluid and the structure at the instant t n ; F n+1 and S n+1 are the status at the instant t n+1 .

①
As we can observe that in such a partitioned procedure there always exists a time lag between the time integration of fluid and structure. Sometimes this sort of coupling approach is also called the staggered method. Due to this time lag of staggered method, the interface conditions (5.1) or (5.2) can not be enforced, hence such a partitioned procedure is described as a loosely coupled solution algorithm. From the point of view of energy [START_REF] Combescure | A numerical scheme to couple subdomains with different time-steps for predominantly linear transient analysis[END_REF], in the basic partitioned method we inject or dissipate numerically energy at the interface, which lead to degradation of numerical stability property and order of accuracy of the used schemes. Although the stability and accuracy of partitioned methods can be improved by means of prediction techniques, their error remains larger than for a monolithic solution procedure [START_REF] Michler | A monolithic approach to fluid-structure interaction[END_REF]. The approach presented subsequently is classified as a monolithic method which follows a strongly coupled procedure for the numerical simulation of fluid-structure interaction problems.

Zero interface energy condition

As presented in [START_REF] Combescure | A numerical scheme to couple subdomains with different time-steps for predominantly linear transient analysis[END_REF], when coupling two physical domains with two different time integrators, one can preserve the numerical stability as well as the minimal order of accuracy in time for the coupled system, as long as the interface energy is ensured to be zero during the numerical simulation. For example, when one wants to couple two different time integrators which are both second order accurate in time, the coupled system will keep the second order of accuracy, if neither energy injection nor energy dissipation occurs at the interface. In contrast, if the zero interface energy condition cannot be ensured, the coupled system will sometimes possess a first order of accuracy or a zero order of accuracy in time, or even a numerically unstable result. Now let's define the increment of interface energy or simply interface energy over the time interval t ∈ [t n , t n+1 ] as

∆W I = t n+1 t n Γ I n s • (-p s I) • v s + n f • (-p f I) • v f dΓ dt (5.4)
where n s • (-p s I) and n f • (-p f I) represent the external forces exerted to the solid and fluid sub-domains in the normal direction, respectively. The product with the velocity gives the power of this external force. Considering the solid and fluid sub-domains as a whole system, we have the total external power occurring at the interface

P I = n s • (-p s I) • v s + n f • (-p f I) • v f (5.5)
then integrating this power along the interface and over the time interval, we obtain ∆W I denoting the external work at the interface. In fact, theoretically this interface energy should be zero, i.e. ∆W I ≡ 0. However, when the numerical methods are used to discretize the space and time, and if the interface conditions are not wisely treated, there will be production or dissipation of algorithmic energy at the interface. Concretely, if ∆W I > 0, we inject energy into the system; if ∆W I < 0, we dissipate energy from the system.

Remark: In this interface energy, only the work due to the normal force has been considered, which corresponds to the slipping condition at the interface with an inviscid fluid.

Consider that the interface

Γ I is discretized into N k elements, i.e. Γ I = N k k Γ k .
In 3-D cases, each element Γ k is assumed to be a planar surface, as shown in Figure 5.3; in 2-D cases, Γ k is described by a section of line; in 1-D cases, it is just a point.

With such an approximation, the normal vector n k is constant inside each element Γ k . Furthermore, for the sake of simplicity, let's assume that the pressure and the velocity fields at the interface are also piece-wise uniform. Now we give the first interface condition used in this work, the equilibrium condition

p f (x k ) = p s (x k ) = p k (5.6)
where x k denotes the position of the interface element Γ k , which may be the geometric center of the element. During the numerical simulation, we will ensure this equilibrium condition in a strong sense, which means that this equation is verified at each instant, i.e.

p f (x n k ) = p s (x n k ) = p n k , p f (x n+1 k ) = p s (x n+1 k ) = p n+1
k and so on.

Then, the spatially discretized interface energy is calculated as

∆W d I = t n+1 t n N k k p k n k • v s (x k ) -v f (x k ) s k dt (5.7)
where n k = n f (x k ) = -n s (x k ) denoting the normal vector of the interface element, and s k is the surface area.

Notice that the interface pressure p k , the normal vector n k , the surface area s k and the velocities are all time-depending variables, which makes the calculation of time integration (5.7) quite difficult to be carried out. To simplify this time integration, some approximation or assumption may be used. For instance, firstly, let's assume that the interface pressure p k is piece-wise constant in time and equals the mean value p k over the time interval t ∈ [t n , t n+1 ], with

p k = 1 ∆t t n+1 t n p k (t) dt (5.8)
With this assumption, the discretized interface energy (5.7) becomes

∆W d I = N k k p k t n+1 t n n k • v s (x k ) -v f (x k ) s k dt (5.9)
Similarly, we can define the mean value of

n k • v s (x k ) -v f (x k ) s k as n k • v s (x k ) -v f (x k ) s k = 1 ∆t t n+1 t n n k • v s (x k ) -v f (x k ) s k dt (5.10)
with which we can rewrite the discretized interface energy as

∆W d I = N k k p k ∆t n k • v s (x k ) -v f (x k ) s k (5.11)
Hence, if we want to ensure the zero interface energy condition, it is sufficient to impose at each surface element that

n k • v s (x k ) -v f (x k ) s k = 0 (5.12)
which is the key point of the coupling method presented in the PhD work. Now we show the different treatments to this condition (5.12) in different cases:

• In geometrically linear cases, the configuration of the structure will not change, hence it is allowed to suppose that the normal vector n k and the surface area s k keep constant as time varies. In such linear cases, (5.12) becomes

n k • v s (x k ) -v f (x k ) s k = 0 (5.13)
hence, we can ensure ∆W d I = 0 by imposing the equality of the mean velocity in the normal direction

n k • v s (x k ) -n k • v f (x k ) = 0 (5.14)
which may be written as

n k • v s (x n+1 k ) + v s (x n k ) 2 -n k • v f (x n+1 k ) + v f (x n k ) 2 = 0 (5.15)
• In geometrically non-linear cases, the shape of the structure will change in time. So the constant geometry assumption is no more valid. As n k and s k are timevarying, the zero interface energy condition ∆W d I = 0 can be ensured by imposing

n k • v s (x k )s k -n k • v f (x k )s k = 0 (5.16)
which may be approximated by

1 2 n n+1 k • v s (x n+1 k )s n+1 k + n n k • v s (x n k )s n k - 1 2 n n+1 k • v f (x n+1 k )s n+1 k + n n k • v f (x n k )s n k = 0 (5.17)
Finally, we obtain the second interface condition for the continuity of velocity in a weak sense, in geometrically linear cases (5.15) as well as in geometrically non-linear cases (5.17). If the two interface conditions are satisfied at the same time, the interface energy can be rigorously ensured to be zero during the numerical simulation, thus the numerical stability property and the order of accuracy in time can be then retained.

The coupled equations to solve

This proposed coupling method is classified as monolithic procedure, because the status of variables in one sub-domain cannot be updated to next time step without considering the interaction effect from the other sub-domain. This feature may retain the numerical stability as well as the order of accuracy in time of used schemes, however in such a monolithic method we have to modify the solution procedures for the solid and the fluid to a certain extent.

Consequence for the fluid

In Chapter 3, we have seen that the semi-discrete system of fluid equations (3.94) will be resolved in time by means of the explicit Runge-Kutta numerical scheme. In a pure fluid problems, at the solid wall boundaries, the velocity of solid wall is often given as an already known condition. By using the method based on resolving a partial Riemann problem at the solid boundary, we can update the state of the volume ω i and the momentum ρ i ω i v i for the fluid particle i, with the help the imposed solid boundary velocity v k . But in a fluid-structure interaction problem, the velocity of the solid boundary cannot be given explicitly for the fluid to advance in time, because the value of the velocity is not imposed but is the solution of the coupled problem.

The partial Riemann problem at the solid boundary offers a relation between the fluid pressure at the solid wall and the velocity of this solid elemental surface. With (3.79) and (3.87), we obtain

p k = i∈D k ω i W ik 2[p i -ρ i c i (v k -v i ) • n k ] (5.18)
where v k = v f (x k ) representing the fluid velocity at the surface element x k . Now let's note v fB,k ≡ v k • n k , then (5.18) can be rewritten as

i∈D k 2ω i ρ i c i W ik v fB,k + p k = i∈D k 2ω i (p i + ρ i c i v i • n k )W ik (5.19)
which describes a relation between p k and v fB,k for any instant.

Because the second-order explicit Runge-Kutta scheme (mid-point version) contains two stages, the solid boundary conditions should be considered at two instant: t n+1/2 and t n+1 . (Assuming that at the instant t n , all the variables are already known.)

Then writing (5.19) in matrix form for all the interface elements at the two instant, we obtain the system of equations for the state of fluid variables at the interface:

(I-F) System of equations at t n+1/2 for fluid

K c f v n+1/2 fB + Λ n+1/2 = g n f (5.20)
with,

K c f (l, k) =      i∈D k 2ω n+1/2 i ρ n+1/2 i c n+1/2 i W n+1/2 ik if l = k 0 if l = k (5.21)
, and

           v n+1/2 fB (k) = v n+1/2 fB,k Λ n+1/2 (k) = p n+1/2 k g n f (k) = i∈D k 2ω n+1/2 i (p n+1/2 i + ρ n+1/2 i c n+1/2 i v n+1/2 i • n n+1/2 k )W n+1/2 ik (5.22)
(II-F) System of equations at t n+1 for fluid

K c f v n+1 fB + Λ n+1 = g n f (5.23) with, K c f (l, k) =      i∈D k 2ω n+1 i ρ n+1 i c n+1 i W n+1 ik if l = k 0 if l = k (5.24)
, and

           v n+1 fB (k) = v n+1 fB,k Λ n+1 (k) = p n+1 k g n f (k) = i∈D k 2ω n+1 i (p n+1 i + ρ n+1 i c n+1 i v n+1 i • n n+1 k )W n+1 ik (5.25)
These two system of equations (5.20) and (5.23) can only give a relation between the interface pressure and the normal fluid velocity at the interface. Hence the equations of the solid sub-domain as well as the interface are needed to resolve completely the FSI problem.

Consequence for the solid

As presented previously, the second-order mid-point version of Runge-Kutta scheme requires two stages, hence let's write the momentum equation (4.31) for the solid subdomain at two instant, t n+1/2 and t n+1 :

M s a n+1/2 s = f n+1/2 ext -f n+1/2 int M s a n+1 s = f n+1 ext -f n+1 int (5.26) Then let's define X n+1/2 ≃ X = (X n+1 +X n )/2
, where X may be the displacement u s , the velocity v s , or the acceleration a s . Combining the Newmark scheme (4.60) with each of the two equations in (5.26), we can obtain two systems of equations for the solid sub-domain, for example, in linear geometry cases [START_REF] Li | Coupling of finite volume and finite element subdomains using different time integrators[END_REF]49]:

(I-S) System of equations at t n+1/2 for linear structure

K c s v n+1/2 s + L n+1/2 p Λ n+1/2 = g n s (5.27) with,                K c s = 1 ∆tγ M s + ∆tβ γ K s g n s = M s ( 1 ∆tγ v n s + 1 -2γ 2γ a n s )- K s u n s + ∆t( 1 2 - β γ )v n s + ∆t 2 ( 1 4 - β 2γ )a n s + f n+1/2 ext,s (5.28)
(II-S) System of equations at t n+1 for linear structure

K c s v n+1 s + L n+1 p Λ n+1 = g n s (5.29) with,                K c s = 1 ∆tγ M s + ∆tβ γ K s g n s = M s ( 1 ∆tγ v n s + 1 -γ γ a n s )- K s u n s + ∆t(1 - β γ )v n s + ∆t 2 ( 1 2 - β γ )a n s + f n+1 ext,s (5.30)

Once v n+1

s is obtained, we can update the displacement and acceleration vectors for the next time step by using

       u n+1 s = ∆tβ γ v n+1 s + u n s + ∆t(1 - β γ )v n s + ∆t 2 ( 1 2 - β γ )a n s a n+1 s = 1 ∆tγ v n+1 s - 1 ∆tγ v n s - 1 -γ γ a n s (5.31)
In (5.27) and (5.29), the matrix L p is a geometry operator which relates the interface pressure field Λ with the fluid-inducing external nodal forces f ext,f . In linear geometry cases, this matrix can be considered to be constant, since the geometry of the structure will not change too much. Usually this geometric non-linearity requires a iterative solution procedure, which makes the FSI problem numerically more difficult to resolve. In this PhD work, for not iterating at each time-step, or even each sub-time-step, we choose to use the explicit Newmark scheme (β = 0 and γ = 0.5) when the deformation of the structure is too large to be neglected. With such an explicit scheme, the displacement field can be calculated exactly by

       u n+1/2 s = u n s + ∆t 2 v n s + ∆t 2 4 a n s u n+1 s = u n s + ∆tv n s + ∆t 2 2 a n s (5.32)
with which the geometry of the structure can be predicted independently of the FSI effects, hence we can obtain the geometry operator L n+1/2 p and L n+1 p . Furthermore, if the structural material property is non-linear, the internal nodal forces can no more be calculated by f int = K s u s . In such a non-linear case, the same explicit Newmark scheme will be applied to obtain firstly the displacement fields u n+1/2 s and u n+1 s , which is used to calculate the strain field and finally the internal nodal forces f n+1/2 int and f n+1 int with the according constitutive modes. Hence, in non-linear cases, the system of equations for the solid sub-domain can be written in the similar form as the linear one, except that (5.28) changes to (NL)

       K c s = 1 ∆tγ M s g n s = M s ( 1 ∆tγ v n s + 1 -2γ 2γ a n s ) -f n+1/2 int + f n+1/2 ext,s (5.33)
and (5.30) is replaced by

(NL)        K c s = 1 ∆tγ M s g n s = M s ( 1 ∆tγ v n s + 1 -γ γ a n s ) -f n+1 int + f n+1 ext,s (5.34)

The coupled system of equations

In the coupling strategy presented in this work, there are two sub-steps for one time step.

For each sub-step, we will solve a system of equations in order to obtain the interface status at the corresponding instant.

• For the instant t = t n+1/2 (I)

       (I-S): K c s v n+1/2 s + L n+1/2 p Λ n+1/2 = g n s (I-F): K c f v n+1/2 fB + Λ n+1/2 = g n f (I-C): L n+1/2 s v n+1/2 s + v n+1/2 fB = 0 (5.35)
• For the instant t = t n+1

(II)

       (II-S): K c s v n+1 s + L n+1 p Λ n+1 = g n s (II-F): K c f v n+1 fB + Λ n+1 = g n f (II-C): L n+1 s v n+1 s + v n+1 fB = -L n s v n s -v n fB (5.36)
where (I-C) and (II-C) represent the velocity condition (5.14) or (5.16) at the interface, which can ensure the zero interface energy condition. L s denotes another geometry operator that links the normal solid velocity at each interface elements to the solid velocity vector v s . In linear geometry cases, L s is constant in time, whereas in non-linear geometry cases, L n+1/2 s and L n+1 s are calculated with u n+1/2 s and u n+1 s , respectively, obtained by using the explicit Newmark scheme at the begin of each sub-step.

As for the solution procedure of this coupled system of equations, we will apply the method used in [49] which is based on the Schur's dual formulation. In (5.35) and (5.36), the unknowns to be calculated are v s , v fB and Λ. Hence, in general the system of equations can be written as

       K c s v s + L p Λ = g n s K c f v fB + Λ = g n f L s v s + v fB = W n (5.37)
where W n = 0 for the instant t n+1/2 , and

W n = -L n s v n s -v n
fB for the instant t n+1 . To resolve such a system of equations, five stages are needed:

• Stage 1: Calculate the free velocity

v f ree s = [K c s ] -1 g n s v f ree fB = [K c f ] -1 g n f (5.38) • Stage 2: Calculate the condensed matrix H H = L s [K c s ] -1 L p + [K c f ] -1 (5.39)
• Stage 3: Calculate the interface pressure field

Λ Λ = [H] -1 {L s v f ree s + v f ree fB -W n } (5.40)
• Stage 4: Calculate the link velocity

v link s = [K c s ] -1 (-L s Λ) v link fB = [K c f ] -1 (-Λ) (5.41)
• Stage 5: Calculate the complete velocity

v s = v f ree s + v link s v fB = v f ree fB + v link fB (5.42)

The coupling algorithm

To offer a survey of the proposed coupling method, we give the coupling procedure in Figure 5.4, and the coupling algorithm is presented as follows:

I-1: The fluid solver receives the already known interface status U n I , and then calculates the fluid status U n+1/2 f ; I-2: The coupler uses U n+1/2 f and the solid status U n s to solve the system of equations (5.35) in order to obtain the interface status U

Chapter 6 Numerical examples

We have proposed one method for the interface conditions to ensure the zero interface energy condition in order to preserve the numerical stability property and the optimal order of accuracy in time of the used numerical schemes in all the sub-domains. So far this energy-conserving coupling method is not yet verified or validated. For verifying a numerical method, we can compare the obtained numerical solution with the according analytical or exact solution, if it exists, or we can also compare with other numerical solutions etc. The verification process is aim at verifying whether the tested numerical method resolves the governing equations right, whereas the validation of a numerical method consists in validating whether the numerical method solves the right problem [START_REF] Roache | Verification and Validation in Computational Science and Engineering[END_REF]. Usually we need to compare the numerical solution with the experimental results in the validation tests. In this work, we did not carry out any experimental validation, but only the comparison with analytical solution and other numerical results. Hence the words verification and validation have no difference in this chapter. This chapter starts with the validation of the energy-conserving coupling method in time integration aspects, with the numerical methods using computational meshes for both structure and fluid sub-domains. Then the numerical method in fluid sub-domain is replaced by the SPH-ALE method, and some 1-D numerical examples are given. Two 2-D numerical test cases are then shown with the presence of free surface flows and large structural deformation. Finally, the result of a 3-D example is given and compared with other coupling methods.

1-D validation test cases

For validating the proposed coupling method, beginning with a simple test case seems to be a wise choice, because the numerical simulation is easy to carry out and usually we can compare the numerical results with the exact solution of this relatively simple test problem. Hence in this section, we will give the results of some 1-D numerical examples as well as the comparison with the analytical solution.

Time integration aspects

Before entering the part of validation of the proposed method for coupling of the SPH-ALE method and the FE method, let's see some numerical examples which are used to validate the idea of conserving the interface energy, by means of numerical methods using computational mesh in both the structure and the fluid sub-domains. In this subsection, we only show some important results and some key points, more attention is given on the analysis of the obtained results. More details can be found in the work of Z. Li [START_REF] Li | Coupling of finite volume and finite element subdomains using different time integrators[END_REF].

x Figure 6.1: Mass-spring system coupled with column of acoustic fluid.

The first example is about the mass-spring system coupled with an acoustic fluid as shown in Figure 6.1. The fluid sub-domain is discretized by the Finite Difference (FD) method with a node-centered mesh (Figure 6.2). The linearized Euler equations are used as the governing equations for the fluid sub-domain in linear gas dynamics problems.

The zero interface energy condition is ensured by imposing

v s = v f (6.1)
x Figure 6.2: Node-centered mesh for the fluid.

or,

v n+1 s + v n s 2 = v n+1 f + v n f 2 (6.2)
where v s denotes the velocity of the mass point, and v f is the fluid velocity at the moving solid boundary.

Remark: In such an energy-conserving coupling method, we impose the velocity continuity condition in a weak sense which ensures rigorously that neither energy injection nor energy dissipation will occur at the fluid structure interface. It is one of the attractive points of the proposed method, since the velocity condition at the interface may be imposed more flexibly, which allows the use of incompatible time steps in different sub-domains [START_REF] Li | Coupling of finite volume and finite element subdomains using different time integrators[END_REF]49].

As for the time integrators in each sub-domain, the second-order Newmark scheme (β = 0.25 and γ = 0.5) is used for the solid system, and two different numerical schemes are used: the Lax-Wendroff (L-W) scheme [START_REF] Hirsch | Numerical computation of internal and external flows[END_REF] which is second-order accurate in time and in space; the 4 th order Runge-Kutta explicit time integrator combining with central difference scheme for the space, which is noted as RK4-CD in the following parts.

Given an initial status of the coupled system Figure 6.3 shows the evolution in time of the displacement u s (t), the velocity v s (t) and the acceleration a s (t) of the mass point, which are calculated with two different numerical schemes for the fluid sub-domain. The numerical results are also compared with the analytical solution of this coupled problem. The procedure for obtaining the analytical solution is given in Appendix A. We can observe that the numerical results are in very good agreement with the analytical solution. Moreover, to verify if the proposed coupling strategy can preserve the order of accuracy, a study of convergence rate is carried out and shown in Figure 6.4 . The error is defined as ERROR

Fluid      ρ f (x, 0) = ρ ref f v f (x, 0) = sin[ ω 1 c ref f (L f -x)] + sin[ ω 2 c ref f (L f -x)] (6.3) Solid        u s (0) = 0 v s (0) = v f (0, 0)
= X num -X ana ∞ X ana ∞ (6.5)
where X = u s , v s or a s and ∞ denotes the maximum norm.

Since the CFL number plays a very important role on dispersion and dissipation errors, we will fix the CFL number while changing the time step. The different time Table 6.1: Different time steps used for the analysis of convergence rate.

Computation N o 1 N o 2 N o 3 N o 4
∆t (s) 0.5 × 10 -5 1 × 10 -5 2 × 10 -5 4 × 10 -5 ∆x (m) 0.25 × 10 -2 0.5 × 10 -2 1 × 10 -2 2 × 10 -2 CFL 0.6564 0.6564 0.6564 0.6564 steps and the corresponding meshes are shown in Table 6.1. The total period of the calculation is 0.02s for the four cases. In Figure 6.4, we can observe that when coupling the 2 nd order Newmark time integrator with 2 nd order Lax-Wendroff scheme or the 4 th order Runge-Kutta method, the numerical results of the coupled system possess a 2 nd order accuracy in time. This can prove that coupling of different time integrators with zero interface energy can preserve the minimal order of accuracy in time of the used numerical schemes in all the sub-domains.

Furthermore, this energy-conserving coupling method can also be used for coupling of finite volume method in ALE settings for fluid and finite element method for structure. Since it is not the main objective of this PhD work, we will not show the results in this chapter. More information can be found in [START_REF] Li | Coupling of finite volume and finite element subdomains using different time integrators[END_REF].

1-D mass-spring system coupled with a column of fluid

So far, we have validated the proposed energy-conserving coupling method by means of some numerical tests using computational mesh for fluid. Now, let's replace the numerical method for the fluid by the SPH-ALE method. The configuration of the first test case is presented in Figure 6.5. The mass-spring system is resolved in time with the Newmark time integrator (β = 0.25 and γ = 0.5). The fluid sub-domain is discretized with the SPH-ALE method in Lagrangian mode, the semi-discrete system of equations is integrated with the 2 nd order Runge-Kutta numerical scheme.

The objective of this linear test case is to verify if the coupling strategy can preserve the order of accuracy in time. For that purpose, a parametric study of the case was performed with a varying time step. The time step was chosen as ∆t = τ, 2τ and 4τ , with τ = 0.25 × 10 -6 s, and the rate of convergence of the numerical solution of the coupled problem was analyzed. Since the used time integrators for the solid and fluid domains are both second order accurate in time, hence if one ensures the zero interface energy condition, one should obtain a second order accuracy in time for the coupled system. To verify numerically the order of accuracy in time, we choose to apply the method used in [START_REF] Michler | A monolithic approach to fluid-structure interaction[END_REF], which determines the observed order of accuracy by calculating

θ = ln X num 4τ -X num 2τ X num 2τ -X num τ / ln(2) (6.6)
where θ denotes the observed order of accuracy in time, X a certain variable calculated by using different time steps (τ , 2τ and 4τ ) and means "L 2 -norm" or "∞-norm". Table 6.2 shows the result of the convergence rate analysis for the displacement of the structure, the density and velocity of the fluid particle which is initially located at x = 0.25L f . One observes that the coupled system possesses a second order of accuracy in time. In the second test case, we replace the mass-spring system by a 1D linear beam shown in Figure 6.7. The initial length L 0 s = 1 m, the initial solid density ρ 0 s = 2700 kg/m 3 , and the initial section area A 0 s = 0.01 m 2 . The Young's modulus E s = 67.5 GPa. The solid beam is discretized into 200 linear truss finite elements, and ∆t = 10 -6 s. A uniform pressure step of 20 MPa is imposed at the time t = 0 s in the fluid cavity, as shown in Figure 6.8. The initial discontinuity is located at the fluid-structure interface. The total period of calculation T = 100∆t. In Figure 6.9 and Figure 6.10, the numerical results obtained with the implicit Newmark scheme (β = 0.25 and γ = 0.5) and the explicit Newmark scheme (β = 0 and γ = 0.5), respectively, are compared with the analytical solution obtained by applying the method used in [START_REF] Liu | Numerical simulation of fluid-structure interaction using modified ghost fluid method and naviers equations[END_REF], which is also given in Appendix B.

One can observe that the comparison is quite good. In both of the two results, the level of the pressure and velocity, and the continuity of these values are correctly calculated, comparing with the analytical solution. With the implicit Newmark scheme for the structure, we got some numerical oscillation, but the simulation is stable. Let us also quote that with the chosen time step and a constant element length the explicit time integrator gives the exact solution in the solid. The induced numerical dissipation is due to the numerical schemes used in the fluid sub-domain, i.e. the use of Riemann solver for the SPH method combined with the Runge-Kutta time integrator. V (m/s)

V f ana V f ana V f num V f num (b)
Figure 6.9: Comparison between the analytical solution and the numerical results using the 1 st initial configuration with the implicit Newmark scheme.

Next, we consider another initial configuration which is shown in Figure 6.11-(a). The initial pressure jump is given within the fluid for x > 0.2L f . The numerical result is also compared with the analytical solution, and a good agreement can be observed. As shown in Figure 6.11, at the beginning the pressure discontinuity created a shock wave propagating to the left and a rarefaction wave propagating to the right; the subfigure (c) represents the theoretical instant for the front of the shock wave to impact on the interface, we observe that due to the numerical dissipation error of the SPH-ALE V (m/s)

V f ana V s ana V f num V s num (b)
Figure 6.10: Comparison between the analytical solution and the numerical results using the 1 st initial configuration with the explicit Newmark scheme.

method, the profile of pressure and velocity are not very sharp at the interface; then this shock wave is divided into an incident wave entering in the structure, and a reflection wave which comes back to the fluid sub-domain, as shown in the sub-figure (d).

Even though some numerical error appear in the results, with the coupling strategy presented in this work, one can correctly calculate the interface status when a shock wave impacts on the fluid structure interface. And it is shown that the coupling method does not influence too much the speeds of the incident and the reflection waves.

1-D propagation of shock wave -non-linear structure

Now, we replace the linear beam by a non-linear behavior beam. The elasto-plastic material property is shown in Figure 6.12 with the Young's modulus E le s = 67.5 GPa, the tangent modulus E lp s = 0.1 × E le s and the yield stress σ 0 = 0.0001 × E le s . In Figure 6.13, we compare the numerical results with the ones obtained by using the linear beam. One can observe that the incident shock wave divides into two waves: a elastic wave and a plastic wave, the former moves faster than the latter Whereas, in the linear beam case, there is only one elastic wave in the structure. Figure 6.12: Elasto-plastic material property.

2-D examples

One of the advantages of using the SPH method in FSI problems is that we can easily track the movement of the fluid-structure interface or the geometric outline of the structure. As such, in this section, we will present two 2-D examples of numerical simulation with the presence of free surface flows interacting with a solid structure under large displacement or deformation.

Notice that in such non-linear cases, the geometry of the structure will vary significantly over time. Hence, the coupling strategy and equations, see (5.33) and (5.34), will be applied for treating the considerable change of structural geometry.

Dam-break flow through an elastic gate

The initial configuration of this test case is shown in Figure 6.14. At the moment t = 0 s, the water is contained in a rigid wall tank which has an elastic rubber gate clamped on its top. The whole system is initially at rest and the water is in hydrostatic equilibrium. The background pressure is 0 Pa, which implies the zero pressure condition at the free surface of the fluid. The geometric dimensions of the system and the discretization parameters are shown in Table 6.3. The elastic gate is discretized into 395 quadrangular elements, with each element possessing 4 nodes, and we use 14000 particles to simulate the fluid domain. As shown in Figure 6.14, a pressure gradient is present in the water along the vertical direction due to the gravity (g = -9.81e z m/s 2 ). This pressure gradient makes the fluid tend to flow out of the water tank by pushing away the elastic gate.

We have carried out a numerical simulation of this FSI problem for a period of 0.4 second. The fluid and solid solvers used the same time-step ∆t = 5 × 10 -6 second. The results of the simulation are shown in Figure 6.16 for every 0.04 second from t = 0 s. Two CPUs of 3.0 GHz are used for this numerical simulation, one CPU for the solid solver, one for the fluid solver. The numerical computation costed 17 h for 80000 timesteps. toci et al. [START_REF] Antoci | Numerical simulation of fluidstructure interfaction by SPH[END_REF] and Rafiee et al. [START_REF] Rafiee | An SPH projection method for simulating fluid-hypoelastic structure interaction[END_REF] have proposed the numerical simulation methods for this sort of FSI problems by using the SPH method for both the fluid and the solid domain. Antoci et al. [START_REF] Antoci | Numerical simulation of fluidstructure interfaction by SPH[END_REF] also compared the numerical solution with the experimental result. Yang et al.

[83] and Fourey et al. [START_REF] Fourey | Violent fluid-structure interaction simulations using a coupled SPH/FEM method[END_REF] have carried out the numerical simulations by coupling the SPH method for the fluid with the FE method for the structure, which is quite similar with the proposed method in this thesis, however, the main differences may be found in two aspects: (1) we use the different methods for calculating fluid pressure at solid wall; (2) their methods are both classified as partitioned coupling procedures which have sometimes the numerical instability problems, because such type of coupling strategy can not ensure the zero interface energy condition. A comparison of the numerical results is shown in Figure 6.17. As one can observe that the horizontal and vertical displacement of this free end of the elastic gate is well calculated with the proposed coupling method, compared with the experimental results. However, it can be also seen that in our numerical result, this elastic gate has been opened a little bit earlier. The reason could be various, since it is a Fluid-Structure Interaction problem, in which any inaccuracy inside any subdomain may cause such disagreement. Nevertheless, we found that the key factor is the material model used for the elastic gate, which means that, the better the numerical material model fits the experimental data in the tensile test, the better numerical result we can obtain in this 2-D numerical test.

Breaking dam on an elastic wall

Next, we investigate a second 2-D test case of fluid-structure interaction problem, of which the initial configuration is shown in Figure 6.18. As one can observe, in a rigid wall container a column of water is initially located at the left side, which is in hydrostatic equilibrium. An elastic wall is clamped in the middle of the bottom of the rigid container wall. The geometric and discretization parameters are given in Table 6.4. The geometry of the system is the same with the one used by Walhorn et al. [START_REF] Walhorn | Fluidstructure coupling within a monolithic model involving free surface flows[END_REF]. They have studied numerically this FSI problem by means of the Level-Set method combined with moving mesh technique. In this thesis, we have taken the same geometric dimensions for the system. Meanwhile, we have chosen the same material property for the elastic wall. The initial solid density ρ 0 s = 2500 kg/m 3 , the Young's modulus E s = 10 6 Pa and the Poisson's ratio ν = 0. As Walhorn et al. [START_REF] Walhorn | Fluidstructure coupling within a monolithic model involving free surface flows[END_REF] did, we applied a linear elasticity model which gives a linear relationship between the Green strain tensor E and the second Piola-Kirchhoff stress tensor S S = C : E

where C is a fourth-order tensor of elastic moduli which is invariant over time. Such a material is called a Saint Venant-Kirchhoff material, or Kirchhoff material for brevity [START_REF] Belytschko | Nonlinear Finite Elements for Continua and Structures[END_REF]. Subsequently, we will compare our numerical result with the one obtained by other authors. In this test case, we also used the same time-step ∆t = 2 × 10 -5 s for the fluid and solid solvers. 50000 time-steps have been carried out for a 1 second physical period. Just as we did in the previous test case, two CPUs were used for the coupling simulation (one for the fluid solver, one for the solid solver). The computation time is 13 h.

In Figure 6.19 one can find the time-history of the horizontal displacement of the upper left corner of the elastic wall, which is compared with the results of Walhorn et al. [START_REF] Walhorn | Fluidstructure coupling within a monolithic model involving free surface flows[END_REF], Idelsohn et al. [START_REF] Idelsohn | Unified Lagrangian formulation for elastic solids and incomressible fluids: Application to fluid-structure interaction problems via the PFEM[END_REF] and Rafiee et al. [START_REF] Rafiee | An SPH projection method for simulating fluid-hypoelastic structure interaction[END_REF]. Figure 6.20 shows the fluid pressure field P f and the structural stress component σ xx of this FSI problem for every 0.05 second.

As the proposed coupling method is supposed to conserve the energy at the fluidstructure interface, i.e. the coupling procedure should not introduce numerically energy dissipation or energy injection at the interface, we verify such feature by giving the mechanical energy balance for each sub-domain in Figure 6.21 and Figure 6.22 as well as an energy balance analysis for the whole coupled system in Figure 6.23.

As one can observe in the energy balance analysis for the structure part (Figure 6.21), the kinetic energy E kin and the deformation potential energy E def start to increase around t = 0.2 s due to the impact from the breaking dam flow. Whereas, from this moment, the gravitational potential energy has slightly decreased, because the fluid flow makes the elastic wall to bend towards to the x direction, which induces a lower position 105 (i)

(j) (k) (l) (m) (n) (o) (p)
d I = n t n+1 t n N k k p k n k • v s (x k ) -v f (x k ) s k dt = n ∆t N k k p k n k • v s (x k ) -v f (x k ) s k ≃ n N k k ∆t 2 p n+1 k n n+1 k • v s (x n+1 k ) -v f (x n+1 k ) s n+1 k + p n k n n k • v s (x n k ) -v f (x n k ) s n k (6.8)
where ∆W d I represents the interface energy E I shown in Figure 6.23. It is obvious that the interface energy is rigorously zero during the numerical simulation, which means that the proposed coupling procedure induces neither energy injection nor energy dissipation at the fluid-structure interface. Therefore, we can conclude that the dissipation of the total energy (E T in Figure 6.23) is purely numerical and entirely introduced by the numerical scheme used in the fluid sub-domain. And more importantly, it turns out that the proposed coupling method is energy-conserving.

Remark: When large displacement or deformation appears for structure, as shown in Figure 6.24, the resolution of equations is a non-linear procedure, since the internal nodal force f n+1 int can no longer be linearly related with the displacement field u n+1 s , in addition, the external nodal force f n+1 ext depends on the new structural geometry. Hence, we linearize this procedure by using the explicit Newmark time integrator. With such an explicit scheme, we can calculate the new displacement field u n+1 s with the physical quantities we have already obtained at the previous time-step. Then this new displacement field will be given to the solid solver as an imposed-displacement condition. The solid solver can give the internal nodal force f n+1 int depending on u n+1 s . Meanwhile, with the new displacement field, we can update the structural geometry for the new moment. Figure 6.24: Fluid-structure interaction problem with large displacement or deformation for the structure.

3-D example -linear structural response

A 3D example is given to show the feasibility for multi-dimensional problems of the proposed coupling approach. Figure 6.25 shows a water filled tank which is made of five rigid walls and one deformable square plate clamped at the four sides. The thickness of the plate is d = 0.01 m, the side of the square is a = 0.1 m, the Young's modulus E s = 100 GPa, the density ρ s = 2500 kg/m 3 , and the Poisson's ratio ν = 0. Here, notice that in our calculation, negative pressure is permitted, since we only want to check the quality of the coupling strategy for a linear problem. This may be physically meaningless, but the main objective of this work is to propose and validate numerically a coupling method, rather than to realize the numerical simulation of some test cases which are more realistic. We impose a uniform pulse pressure loading P (t) on the external surface of the deformable plate, which varies in time as presented in Figure 6.27 with p 0 = 0 Pa being equal to the initial pressure of fluid, and p 1 = 10 7 Pa, τ = 0.04 ms. The central point of the plate will have the maximum displacement amplitude in X direction. Figure 6.28 compares the time history of this displacement obtained with the proposed coupling method and the same result obtained with the Europlexus software. Europlexus [START_REF] Europlexus | User's manual[END_REF] is a research finite element explicit code developed by JRC ISPRA and CEA (French Atomic Commission), which uses the usual ALE method for fluid-structure interaction. The total simulation time is 1 ms with 1000 time steps of 1 microsecond. As time evolves, numerical dissipation can be observed for both simulations, but it is higher for the presented method. One will now concentrate on the explanation of this dissipation. Is it due to the presented coupling scheme or to the numerical scheme used for the time integration in each sub-domain? It is known that the mean acceleration Newmark scheme chosen for the solid sub-domain is not dissipative. We have then computed the energies evolution in the whole system and separated the solid, fluid and interface ones. The following equations give the incremental energy variations for the fluid-structure system

                         E S = 1 2 v ⊤ s M s v s + 1 2 u ⊤ s K s u s E F = i 1 2 m i |v f (x i )| 2 + i m i e i E I = n k s k ∆t p n+1 k + p n k 2 n k • v s (x n+1 k ) + v s (x n k ) 2 - v f (x n+1 k ) + v f (x n k ) 2 E T = E S + E F + E I (6.9
) where E S and E F denote the total energy of solid and fluid sub-domain, respectively. They are both the sum of kinetic energy and internal energy. e i is the internal energy per unit mass for the fluid particle i, and calculated by the equation [START_REF] Fourey | Développement d'une méthode de couplage fluide-structure SPH-Éléments Finis en vue de son application à l'hydrodynamique navale[END_REF]:

e i = [p i + (c ref f ) 2 (ρ ref f -ρ i )]/[ρ i (ζ -1)
]. E I is the interface energy and E T is the total energy of the whole coupled system.

-0. The energy balance is presented in Figure 6.29. One can observe that the interface energy is strictly zero during the numerical simulation. The total energy of coupled system increases at first because of the external work of the imposed pressure loading. Then, a decrease of the total energy can be observed clearly, as well as the exchange of energy between the structure and the fluid. It is clear from this analysis that the observed dissipation is due to the used Runge-Kutta time integrator combined with the MUSCL scheme in the fluid. .29: Energy balance of the whole coupled system ("E S ": total energy of structure, "E F ": total energy of fluid, "E I ": interface energy, "E T ": total energy of the whole coupled system). At last, we continue the calculation until 0.004 s with 4000 time steps in order to carry out a spectral analysis of the numerical result of displacement of the center of the plate, obtained by the proposed coupling method. The resulting frequency spectrum is presented in Figure 6.30 At the mean time, we also carry out a modal analysis by using the Cast3M software [START_REF]User's manual[END_REF] on the same problem: analysis of vibro-acoustic modes for coupled fluid-structure system. The whole fluid-structure domain is modeled with 3D 8 nodes brick finite elements (with solid material for the plate and acoustic model for the fluid part). The method used in this formulation is described in [START_REF] Combescure | Fluid-structure interaction: A general method used in the ceasemt computer programs[END_REF]. Table 6.5 compares the first 8 eigenfrequencies of the coupled system to the first 3 eigenfrequencies "measured" with the computed time history of displacement of the center of the plate. Let us also quote that all the low frequency modes which are not present in the spectral analysis of the mid point of the plate correspond to the modes for which the displacement of this node is zero. One observes a very good comparison of the two approaches.

Chapter 7

Conclusions and perspectives

Conclusions

In this PhD thesis, an energy-conserving coupling strategy has been proposed for transient fluid-structure interaction. This coupling method is validated by several 1-D test cases and a 3-D example is also given to show the feasibility for multi-dimensional cases. Through these numerical experiments, we observe that the proposed method can ensure the numerical stability of the coupling simulation, and it can preserve the minimal order of temporal accuracy of the time integrators used in all the sub-domains.

Different numerical methods are used for the two physical domains: the structure is simulated by the finite element method and the fluid sub-domain is discretized into a set of particles by the SPH-ALE method, hence it is easy to track the interface. A velocity constraint condition is imposed in a weak sense at the interface in the normal direction, which can ensure the zero interface condition.

Being different with a staggered coupling procedure, the proposed coupling method is classified as a monolithic approach, in which there is no time lag between the time integration of the two sub-domains. Numerically, neither energy injection nor energy dissipation will occur at the fluid-structure interface.

However, like every monolithic approach, an inconvenient of the proposed coupling method is that some necessary modification of the existing simulation codes will be required. Secondly, the coupling method is applicable for the specific time integrators.

Perspectives

This work includes the theoretical demonstration of the coupling method as well as some important validation test cases. In the future, it will be of great interest to consider the physical viscosity of the fluid. In this thesis, the fluid is assumed to be inviscid, which is described by the Euler equations. Considering the viscous effect, one will resolve the Navier-Stokes equations instead of the Euler equations for the fluid sub-domain. At the interface, the equilibrium and the velocity continuity should be considered not only in the normal direction but also in the tangential direction. As for the numerical simulation aspect, one will resolve some extra coupling equations at the interface in order to take into account the tangential velocity constraint condition.

Secondly, since the velocity constraint condition is imposed in a weak sense, it will be also be interesting to use incompatible time steps to make the numerical simulation more efficient. For example, in rapid dynamic problems, when an explicit time integrator is applied, a small time step is often used for the reason of numerical stability. In such cases, the time step for the solid may be smaller than the one used for the fluid. Hence, one can use different time steps in the numerical simulation of fluid-structure interaction.

since we have T ′′ (t) = -λT (t), hence

-λm s X(L f )T (t) + k s X(L f )T (t) + ρ ref f (c ref f ) 2 A s X ′ (L f )T (t) = 0 (A.14)
Dividing (A.14) by T (t) gives

-λm s X(L f ) + k s X(L f ) + ρ ref f (c ref f ) 2 A s X ′ (L f ) = 0 (A.15)
or,

X ′ (L f ) = m s λ -k s ρ ref f (c ref f ) 2 A s X(L f ) (A.16)
Using the general solution (A.12) in (A. [START_REF] De Leffe | SPH modeling of shallow-water coastal flows[END_REF] gives

C 2 √ λ c ref f cos( √ λ c ref f L f ) = m s λ -k s ρ ref f (c ref f ) 2 A s C 2 sin( √ λ c ref f L f ) (A.17) hence, tan( √ λ c ref f L f ) = √ λρ ref f c ref f A s m s λ -k s (A.18)
By resolving numerically this equation, we can obtain the eigenfrequencies of the coupled system ω n = λ n (A. [START_REF] Dubois | Partial riemann problem, boundary conditions, and gas dynamics. Absorbing Boundaries and Layers, Domain Decomposition Methods: Applications to Large Scale Computations[END_REF] where n denotes the number of the mode. The according eigenfunction write

X n (x) = C n 2 sin( √ λ n c ref f x) (A.20)
Notice that the boundary condition at x = L f is not homogeneous, hence the eigenfunctions are not orthogonal, which means that the inner product of two eigenfunctions X m and X n is not zero when m = n < X m , X n > = L f 0 X m (x)X n (x) dx = 0 (m = n) (A.21)

Therefore, as did in [START_REF] Darmawijoyo | On the weakly damped vibrations of a string attached to a spring-mass-dashpot system[END_REF], we will define a novel inner product under this nonhomogeneous boundary condition. Taking two eigenfunctions X m and X n with two according eigenfrequencies ω m = √ λ m and

ω n = √ λ n          -X ′′ m (x) = λ m (c ref f ) 2 X m (x) -X ′′ n (x) = λ n (c ref f ) 2 X n (x) (A.22)
and we have the identity

-X ′′ m (x)X n (x) + X m (x)X ′′ n (x) = [-X ′ m (x)X n (x) + X m (x)X ′ n (x)] ′ (A.
λ m -λ n (c ref f ) 2 L f 0 X m (x)X n (x) dx = -X ′ m (L f )X n (L f ) + X m (L f )X ′ n (L f ) (A.25)
Substituting (A.16) in (A. [START_REF] Farhat | Robust and provably second-order explicit-explicit and implicit-explicit staggered timeintegrators for highly non-linear compressible fluid-structure interaction problems[END_REF] gives

λ m -λ n (c ref f ) 2 L f 0 X m (x)X n (x) dx = - m s (λ m -λ n ) ρ ref f (c ref f ) 2 A s X m (L f )X n (L f ) (A.26) 123
where X m (L f )X n (L f ) can be expressed by

X m (L f )X n (L f ) = L f 0 δ(x -L f )X m (x)X n (x) dx (A.
λ m -λ n (c ref f ) 2 L f 0 X m (x)X n (x) dx = - λ m -λ n (c ref f ) 2 L f 0 m s ρ ref f A s δ(x -L f )X m (x)X n (x) dx (A.28) or, λ m -λ n (c ref f ) 2 L f 0 1 + m s ρ ref f A s δ(x -L f ) X m (x)X n (x) dx = 0 (A.29)
which implies that when m = n hence,

L f 0 1 + m s ρ ref f A s δ(x -L f ) X m (x)X n (x) dx = 0 (A.
L f 0 1 + m s ρ ref f A s δ(x -L f ) sin( √ λ m c ref f x) ∞ n=1 B n sin( √ λ n c ref f x) dx = L f 0 1 + m s ρ ref f A s δ(x -L f ) sin( √ λ m c ref f x)φ(x) dx (A.41)
As presented previously, the inner product between two eigenfunctions equals zero when m = n, hence

B m L f 0 1 + m s ρ ref f A s δ(x -L f ) sin 2 ( √ λ m c ref f x) dx = L f 0 1 + m s ρ ref f A s δ(x -L f ) sin( √ λ m c ref f x)φ(x) dx (A.42)
Finally the coefficient B m can be calculated by

B m = L f 0 1 + m s ρ ref f A s δ(x -L f ) sin( √ λ m c ref f x)φ(x) dx L f 0 1 + m s ρ ref f A s δ(x -L f ) sin 2 ( √ λ m c ref f x) dx (A.43)

Appendix B Analytical solution of the 1-D shock wave propagation problem

In this case, we consider a 1-D elastic solid bar connected with a column of water. The behavior of the elastic bar is assumed to be linear, which implies that the deformation of structure is small. Hence, the solid equation writes

ρ s ∂ 2 u s ∂t 2 -E s ∂ 2 u s ∂x 2 = 0 (B.1)
where u s = u s (x, t) denotes the displacement field in the solid sub-domain. ρ s and E s are two constants representing the solid density and the Young's modulus, respectively. The speed of sound is defined as c According to the Generalized Riemann Invariants condition [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics[END_REF], we obtain the two relations across the corresponding waves dσ s + ρ s c s dv s = 0 across the wave associated with λ 1 s = +c s dσ s -ρ s c s dv s = 0 across the wave associated with λ 2 s = -c s (B.10) one of which will be used subsequently to find the interface status in the fluid-structure interaction problem.

As for the fluid part, we have already discussed previously in Section 3. The fluid governing equation writes

∂V f ∂t + A f ∂V f ∂x = 0 (B.11)
with,

V f = ρ f v f and A f = v f ρ f c 2 f /ρ f v f (B.12)
129 which possesses two eigenvalues between which we will chose the appropriate one for the fluid-structure interaction problem.

   λ 1 f = v f + c f λ 2 f = v f -c f (B.
Let's now consider a wave propagation problem across the fluid-structure interface. Suppose that the discontinuity is initially located at the interface. One should resolve a Riemann problem presented in for the part of structure, which corresponds to the relation dσ s -ρ s c s dv s = 0. Since this is a linear elastic wave, we can easily relate the status in the star region with the one in the left region

(σ * s -σ L s ) -ρ s c s (v * s -v L s ) = 0 (B.15)
However, for the fluid part, it is more complicate, because there exists two possibilities for the type of wave:

• If the wave associated with λ 1 f = v f + c f is a shock wave The Rankine-Hugoniot condition gives

   S R f (ρ R f -ρ * f ) = ρ R f v R f -ρ * f v * f S R f (ρ R f v R f -ρ * f v * f ) = ρ R f (v R f ) 2 + p R f -ρ * f (v * f ) 2 + p * f (B.16)
where S R f denotes the speed of wave associated with λ 1 f = v f + c f . And we can obtain

v * f -v R f = (p * f -p R f )( 1 ρ R f - 1 ρ * f ) 1/2
(B.17)

• If the wave associated with λ 1 f = v f + c f is a rarefaction wave We carry out an integration of the corresponding Generalized Riemann Invariants condition across this wave

ρ * f ρ R f c f ρ f dρ f = v * f v R f dv f (B.18)
With the Tait equation (2.44), the speed of sound c f can be expressed as

c f = ∂p f ∂ρ f 1/2 = Bζ (ρ ref f ) ζ ρ ζ-1 f 1/2 (B.19)
with which (B.18) can be rewritten as

Bζ (ρ ref f ) ζ 1/2 ρ * f ρ R f ρ ζ-3 2 f dρ f = v * f -v R f (B.20)
Finally, we obtain the relation across the rarefaction wave

v * f -v R f = Bζ (ρ ref f ) ζ 1/2 2(ρ R f ) ζ-1 2 ζ -1   p * f + B p R f + B ζ-1 2ζ -1   (B.21)
As one can observe that no matter which type of wave, the relation (B.17) and (B.21) are both complicate, and it is hard to resolve the status in the star region with them. For the sake of simplicity for the resolution, we use the method presented by Liu et al. [START_REF] Liu | Numerical simulation of fluid-structure interaction using modified ghost fluid method and naviers equations[END_REF]. We carry out the Taylor development of the right term in (B.17) and (B.21) with respect to p * f , and we only keep the first order term which turn out to be the same for the kinds of wave

v * f -v R f ≃ 1 ρ R f c R f (p * f -p R f ) (B.22)
Gathering (B.15) and (B.22) as well as the continuity condition at the interface, we obtain

                 (σ * s -σ L s ) = ρ s c s (v * s -v L s ) v * f -v R f ≃ 1 ρ R f c R f (p * f -p R f ) v * s = v * f σ * s = -p * f (B.23)
Resolving this system of equations, we can obtain the status of variables in the star region. Then we can determine the type of wave moving into the fluid sub-domain. The subsequent procedures are the same for a classic Riemann problem, and reader might find more details in the book of E.F. TORO [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics[END_REF].

1. 1 5 . 2

 152 Accidents caused by Fluid-Structure Interaction. . . . . . . . . . . . . . 1.2 A typical Pelton turbine. . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 The blades of Pelton turbine. . . . . . . . . . . . . . . . . . . . . . . . 2.1 The mappings between the initial configuration Ω 0 and the current configuration Ω. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Kernel approximation with the smoothing kernel function W (r, h). . . . 3.2 The Riemann problem between the fluid particle i and the particle j. . . 3.3 Structure of the solution of the Riemann problem in the ξ -t plane for the time-dependent Euler equations -three wave families associated with the eigenvalues v ξ -c, v ξ and v ξ + c. . . . . . . . . . . . . . . . . 3.4 Spatial distribution of solution for first order and second order Godunovtype schemes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Boundary conditions treatment for hyperbolic fluid equations. . . . . . 3.6 The truncation of the support domain of the kernel function for the fluid particle i by the solid wall boundary. . . . . . . . . . . . . . . . . . . . 3.7 Partial Riemann problem at the solid wall boundary. . . . . . . . . . . . 3.8 Solving a partial Riemann problem at the solid wall boundary with a imposed value of the density ρ * . . . . . . . . . . . . . . . . . . . . . . 3.9 Advancing in time from the instant t n to t n+1 by using time integrator. . 4.1 The solid domain Ω s discretized by a Lagrangian mesh under two types of boundary conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1 The normal vectors for each sub-domain at the fluid structure interface. ix The procedure of basic partitioned coupling method. F n and S n denote the status of the fluid and the structure at the instant t n ; F n+1 and S n+1 are the status at the instant t n+1 . . . . . . . . . . . . . . . . . . . . . . 5.3 Spatial discretization of the interface Γ I into several elemental surfaces Γ k . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.4 The coupling procedure for fluid-structure interaction. . . . . . . . . . . 6.1 Mass-spring system coupled with column of acoustic fluid. . . . . . . . 6.2 Node-centered mesh for the fluid. . . . . . . . . . . . . . . . . . . . . 6.3 Numerical results with different schemes and analytical solution for the mass-spring system coupled with an acoustic fluid. . . . . . . . . . . . 6.4 Convergence rates of the solid displacement, velocity and acceleration using different schemes for the fluid: (a) L-W; (b) RK4-CD. . . . . . . 6.5 Linear test case -mass-spring system coupled with a column of water. . 6.6 Coupling result of the linear test case -evolution in time of the displacement u s , the velocity v s and the acceleration a s of the mass point (∆t = 10 -6 s). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.7 1D linear beam coupled with a column of water -propagation of shock wave across the fluid-structure interface. . . . . . . . . . . . . . . . . . 6.8 1 st initial configuration for the test case with 1D linear beam -(a) initial profile of fluid pressure (p f ) and solid stress (-σ s ); (b) initial profile of fluid and solid velocities. . . . . . . . . . . . . . . . . . . . . . . . . .

( 2 Ω

 2 nd initial configuration, Implicit Newmark scheme) -(a) t = 0 s; (b) t = 5 × 10 -5 s; (c) t = 1.35 × 10 -4 s; (d) t = 2.35 × 10 -4 s. . . . . . . 6.12 Elasto-plastic material property. . . . . . . . . . . . . . . . . . . . . . 6.13 Coupling results of the 1D beam test case (2 nd initial configuration, t = 2.35 × 10 -4 s) -(a) non-linear beam; (b) linear beam. . . . . . . . . . . x 6.14 Initial configuration of the test case "dam-break flow through an elastic gate". . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.15 The used Mooney-Rivlin material model for non-linear elasticity rubber. 6.16 Numerical simulation results of the test case "dam-break flow through an elastic gate" -fluid pressure field P f and the stress component σ zz in the structure for every 0.04 second. . . . . . . . . . . . . . . . . . . . . 6.16 Numerical simulation results of the test case "dam-break flow through an elastic gate" -fluid pressure field P f and the stress component σ zz in the structure for every 0.04 second. . . . . . . . . . . . . . . . . . . . . 6.17 Comparison of the displacement of the free end of the elastic plate -(a) horizontal displacement; (b) vertical displacement. . . . . . . . . . . . 6.18 Initial configuration of the test case "breaking dam on an elastic wall". . 6.19 Comparison of the upper left corner displacement of the elastic wall between numerical results for the test case "breaking dam on an elastic wall". . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.20 Numerical simulation results of the test case "breaking dam on an elastic wall" -fluid pressure field P f and the stress component σ xx in the structure.107 6.20 Numerical simulation results of the test case "breaking dam on an elastic wall" -fluid pressure field P f and the stress component σ xx in the structure.108 6.21 Mechanical energy balance for the solid sub-domain. . . . . . . . . . . 6.22 Mechanical energy balance for the fluid sub-domain. . . . . . . . . . . 6.23 Mechanical energy balance for the whole coupled system. . . . . . . . . 6.24 Fluid-structure interaction problem with large displacement or deformation for the structure. . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.25 Configuration of the 3D test case with linear structural response. . . . . 6.26 The finite element mesh and the fluid particles. . . . . . . . . . . . . . 6.27 Time varing pressure loading on the external surface of the structure. . . 6.28 Evolution in time of the central point's displacement in X. . . . . . . . 6.29 Energy balance of the whole coupled system ("E S ": total energy of structure, "E F ": total energy of fluid, "E I ": interface energy, "E T ": total energy of the whole coupled system). . . . . . . . . . . . . . . . . . . 6.30 Frequency spectrum of the displacement obtained by using the presented coupling method. . . . . . . . . . . . . . . . . . . . . . . . . . xi A.1 Coordinates system for the problem of mass-spring system coupled with a column of acoustic fluid. . . . . . . . . . . . . . . . . . . . . . . . . B.1 Two characteristic lines in the x -t plane for the solid problem. . . . . B.2 Two non-linear waves in the x -t plane for the fluid problem. . . . . . B.3 Two waves in the x-t plane for fluid-structure interaction problem with the initial discontinuity at the interface. . . . . . . . . . . . . . . . . . xii Nomenclature Symbols Current configuration of a material domain Ω 0 Initial configuration of a material domain χ Reference coordinates Φ Mapping between the initial configuration and the current configuration X Material or Lagrangian coordinates x Spatial or Eulerian coordinates Operators ∇ x Left gradient operator with respect to the spatial coordinates

  (a) Collapse of Tacoma Narrows Bridge (1940). (b) Fracture of water pipe.

Figure 1 . 1 :

 11 Figure 1.1: Accidents caused by Fluid-Structure Interaction.

Figure 1 . 2 :

 12 Figure 1.2: A typical Pelton turbine.Figure 1.3: The blades of Pelton turbine.

Figure 1 . 3 :

 13 Figure 1.2: A typical Pelton turbine.Figure 1.3: The blades of Pelton turbine.

Figure 2 . 1 :

 21 Figure 2.1: The mappings between the initial configuration Ω 0 and the current configuration Ω.

( 2 .

 2 [START_REF] Leduc | Étude physique et numérique de l'écoulement dans un système d'injection de turbine Pelton[END_REF] where ρ f denotes the fluid density, v f the fluid velocity and σ f the tensor of stress in the fluid domain. (2.41) are the conservation equations of the fluid mass and linear momentum within a volume Ω moving with an arbitrary velocity v.

and ζ ≃ 7 ,

 7 with ρ ref f being the reference density and c ref f the reference speed of sound for the fluid.

Figure 3 . 1 :

 31 Figure 3.1: Kernel approximation with the smoothing kernel function W (r, h).

Figure 3 . 3 :

 33 Figure 3.3: Structure of the solution of the Riemann problem in the ξ -t plane for the time-dependent Euler equations -three wave families associated with the eigenvalues v ξ -c, v ξ and v ξ + c.

Figure 3 . 4 :

 34 Figure 3.4: Spatial distribution of solution for first order and second order Godunovtype schemes.

Figure 3 . 5 :

 35 Figure 3.5: Boundary conditions treatment for hyperbolic fluid equations.

Figure 3 . 6 :

 36 Figure 3.6: The truncation of the support domain of the kernel function for the fluid particle i by the solid wall boundary.

( 3 . 70 )

 370 Combining (3.69) and (3.70) gives the semi-discrete equations of (3.35) for the fluid particles near to the solid wall boundary

Figure 3 . 7 :

 37 Figure 3.7: Partial Riemann problem at the solid wall boundary.

. 8 .Figure 3 . 8 :

 838 Figure 3.8: Solving a partial Riemann problem at the solid wall boundary with a imposed value of the density ρ * .

Figure 4 . 1 :

 41 Figure 4.1: The solid domain Ω s discretized by a Lagrangian mesh under two types of boundary conditions.

•

  Virtual internal power δP int δP int = Ωs δD : σ dΩ (4.13) • Virtual external power δP ext δP ext = Ωs δv • ρb dΩ + Γt δv • t b dΓ (4.14) • Virtual kinetic power δP kin δP kin = Ωs δv • ρ v dΩ (4.15) Hence, we can write (4.8) in a concise form δP tot = δP int -δP ext + δP kin = 0 ∀δv ∈ U 0 (4.16)

M 11 IM

 11 and external nodal forces for the node I, respectively, and M IJ is one of the elements in the so called mass matrix M  power principle (4.16) becomesδv ⊤ I M IJ vJ + f int If ext I = 0 (4.28)Because the virtual velocity can be chosen arbitrarily, one can obtain the momentum equation for the node I M IJ vJ + f int Isystem of equations for all the nodes writesM vs + f intf ext = 0(4.30) or, Ma s + f intf ext = 0 M 12 I . . . M 1I I . . . M 1n I M 21 I M 22 I . . . M 2I I . . . M 2n I I1 I M I2 I . . . M II I . . . M In I

  35) is rewritten in matrix form {σ} = [C] {ε} (4.37) with, {σ} = σ xx σ yy σ zz σ xy σ yz σ zx ⊤ (4.38) and, {ε} = ε xx ε yy ε zz 2ε xy 2ε yz 2ε zx ⊤ (4.39)

44 )

 44 Using this matrix B I , we can express the strain tensor {ε} as {ε} = B I u I (4.45)

Figure 5 . 1 :

 51 Figure 5.1: The normal vectors for each sub-domain at the fluid structure interface.

  Step one: the fluid gives the pressure at the instant t n to the structure;②Step two: the structural system advances in time to the instant t n+1 under the given fluid pressure load; ③ Step three: the fluid uses this new position of the interface to update the fluid mesh to the instant t n+1 ; ④ Step four: the fluid system advances in time to t n+1 and obtains a new pressure load for the next step. This procedure is shown in the Figure 5.2.

Figure 5 . 2 :

 52 Figure5.2: The procedure of basic partitioned coupling method. F n and S n denote the status of the fluid and the structure at the instant t n ; F n+1 and S n+1 are the status at the instant t n+1 .

Figure 5 . 3 :

 53 Figure 5.3: Spatial discretization of the interface Γ I into several elemental surfaces Γ k .

2 :;II- 3 :Figure 5 . 4 :

 2354 Figure 5.4: The coupling procedure for fluid-structure interaction.

  /m 3 and c ref f = 328.2 m/s, and ω 1 = 341.6 rad/s and ω 2 = 1172 rad/s. The mass m s = 0.8 kg and the spring stiffness k s = 8000 N/m. The length of tube L f = 1 m.

0 0 .Figure 6 . 3 :

 063 Figure 6.3: Numerical results with different schemes and analytical solution for the mass-spring system coupled with an acoustic fluid.

Figure 6 . 4 :

 64 Figure 6.4: Convergence rates of the solid displacement, velocity and acceleration using different schemes for the fluid: (a) L-W; (b) RK4-CD.

Figure 6 . 5 :

 65 Figure 6.5: Linear test case -mass-spring system coupled with a column of water.

2 )Figure 6 . 6 :

 266 Figure 6.6: Coupling result of the linear test case -evolution in time of the displacement u s , the velocity v s and the acceleration a s of the mass point (∆t = 10 -6 s).

  x

Figure 6 . 7 :

 67 Figure 6.7: 1D linear beam coupled with a column of water -propagation of shock wave across the fluid-structure interface.

Figure 6 . 8 :

 68 Figure 6.8: 1 st initial configuration for the test case with 1D linear beam -(a) initial profile of fluid pressure (p f ) and solid stress (-σ s ); (b) initial profile of fluid and solid velocities.

Figure 6 . 11 :

 611 Figure 6.11: Comparison between the numerical result and the analytical solution (2 nd initial configuration, Implicit Newmark scheme) -(a) t = 0 s; (b) t = 5 × 10 -5 s; (c) t = 1.35 × 10 -4 s; (d) t = 2.35 × 10 -4 s.

Figure 6 . 13 :

 613 Figure 6.13: Coupling results of the 1D beam test case (2 nd initial configuration, t = 2.35 × 10 -4 s) -(a) non-linear beam; (b) linear beam.

Figure 6 . 14 :

 614 Figure 6.14: Initial configuration of the test case "dam-break flow through an elastic gate".

Figure 6 . 15 :

 615 Figure 6.15: The used Mooney-Rivlin material model for non-linear elasticity rubber.

Figure 6 . 16 :Figure 6 . 16 :Figure 6 . 17 :

 616616617 Figure 6.16: Numerical simulation results of the test case "dam-break flow through an elastic gate" -fluid pressure field P f and the stress component σ zz in the structure for every 0.04 second.

Figure 6 . 18 :

 618 Figure 6.18: Initial configuration of the test case "breaking dam on an elastic wall".

Figure 6 . 20 :

 620 Figure 6.20: Numerical simulation results of the test case "breaking dam on an elastic wall" -fluid pressure field P f and the stress component σ xx in the structure.

Figure 6 . 20 :Figure 6 . 21 :Figure 6 .Figure 6 .

 62062166 Figure 6.20: Numerical simulation results of the test case "breaking dam on an elastic wall" -fluid pressure field P f and the stress component σ xx in the structure.

∆W

  3. A structural mesh is used to discretize the plate into 20(a) × 20(a) × 10(d) 8 nodes brick finite elements. The length of the tank L f = 0.2 m, and there are 20 × 20 × 40 fluid particles. Initially, the system is at rest, with the initial fluid pressure equaling the reference pressure, p f (x, t = 0) = p ref f = 0 Pa.

Figure 6 . 25 :

 625 Figure 6.25: Configuration of the 3D test case with linear structural response.

Figure 6 .

 6 Figure 6.26: The finite element mesh and the fluid particles.

Figure 6 . 27 :

 627 Figure 6.27: Time varing pressure loading on the external surface of the structure.

Figure 6 . 28 :

 628 Figure 6.28: Evolution in time of the central point's displacement in X.

Figure 6

 6 Figure 6.29: Energy balance of the whole coupled system ("E S ": total energy of structure, "E F ": total energy of fluid, "E I ": interface energy, "E T ": total energy of the whole coupled system).

Figure 6 . 30 :

 630 Figure 6.30: Frequency spectrum of the displacement obtained by using the presented coupling method.

  23)Then integrating this relation (A.23) over the domainx ∈ [0, L f ] ′′ m (x)X n (x) + X m (x)X ′′ n (x)] dx = [-X ′ m (x)X n (x) + X m (x)X ′ n (x)]L f 0 (A.24) With (A.22) and the boundary condition at x = 0 we have

  [START_REF] Fish | A first course in finite elements[END_REF] with δ(x) being the Delta function.Then (A.26) becomes

30 )

 30 Hence we will define the novel inner product of eigenfunctions as (A.30). With the used initial conditions (6.3), we can obtain that function T (t) in (A.8). The general solution writesT (t) = A 1 cos( √ λt) + A 2 sin( √ λt) (A.32)where A 1 and A 2 are two coefficients to be determined.

  real values, hence, (B.7) is an hyperbolic equation, and can be resolved by using the characteristic method. The two corresponding characteristic lines are shown in Figure B.1.

Figure B. 1 :

 1 Figure B.1: Two characteristic lines in the x -t plane for the solid problem. Remark: Notice that we are going to resolve a linear wave propagation problem, since the eigenvalues are constant. The values of λ 1 s and λ 2 s represent the speed of waves propagating to the right and the left, respectively.

  13) corresponding to the two non-linear waves shown in Figure B.2.

Figure B. 2 :ρ f c f dv f = 0 across the wave associated with λ 1 f = v f + c f dρ f + ρ f c f dv f = 0

 210 Figure B.2: Two non-linear waves in the x -t plane for the fluid problem.

  Figure B.3.

Figure B. 3 :

 3 Figure B.3: Two waves in the x -t plane for fluid-structure interaction problem with the initial discontinuity at the interface.

  s

  δD xx δD yy δD zz 2δD xy 2δD yz 2δD zx

	with,
			⊤
				(4.41)
	Then the virtual internal power δP int becomes
	δP int =	Ωs	{δD} ⊤ {σ} dΩ	(4.42)
	Now let's construct a matrix relating {δD} with the nodal virtual velocity δv I
	{δD} = B I δv I	(4.43)

[START_REF] Johnson | SPH for high velocity impact computations[END_REF]

where E denotes the Young's modulus, and ν is the Poisson's ratio.

Similarly, the virtual strain rate tensor δD in (4.13) can also be expressed with Voigt notation method {δD} =

Table 6 .

 6 2: The observed order of accuracy in time of the coupling result.

	Variables (X)	θ ∞	θ L 2
	u s	2.0483 2.0405
	ρ f (25%)	1.9920 1.9989
	v f (25%)	1.9939 1.9999
	6.1.3 1-D propagation of shock wave -linear structure

Table 6 . 3 :

 63 Dimensions of the system and discretization parameters.

	Water width (A)	0.1 m
	Water depth (H)	0.14 m
	Gate width (s)	0.005 m
	Gate height (L)	0.079 m
	Fluid particle size	0.001 m
	Rigid wall spacing	0.001 m
	Elastic gate mesh spacing	0.001 m

Table 6 .

 6 4: Geometric and discretization parameters.

	Water width (L)	0.146 m
	Water height	2L
	Initial distance	L
	Width of rigid container	4L
	Elastic wall width (b)	0.012 m
	Elastic wall height (a)	20b/3
	Fluid particle size	0.002 m
	Rigid wall spacing	0.002 m
	Elastic wall mesh spacing	0.002 m

Table 6 . 5 :

 65 Comparison of the eigenfrequencies of the coupled system between using the present coupling method and Cast3M.

	F eig (kHz)	f 1	f 2	f 3	f 4	f 5	f 6	f 7
	present method 3.364 6.718	×	×	×	×	8.847
	Cast3M	3.461 6.627 7.328 7.328 8.132 8.132 8.992
	Difference	2.8 % 1.4 %	×	×	×	×	1.6 %

  2 s = E s /ρ s . If we note

	Furthermore, from (B.2) we can also obtain
		      	1 E s	∂v s ∂x ∂σ s ∂t	= =	∂ 2 u s ∂t∂x ∂x∂t ∂ 2 u s	(B.4)
	where the equality between the right terms yields
		∂σ s ∂t	-E s	∂v s ∂x	= 0	(B.5)
	Regrouping (B.3) and (B.5), we obtain the system of equation for the solid sub-
	domain	      	∂v s ∂t ∂σ s ∂t	--E s 1 ρ s ∂σ s ∂x ∂x ∂v s	= 0 = 0	(B.6)
			    	∂u s ∂t ∂x σ s = E s v s = ∂u s	(B.2)
	with v ∂v s ∂t	-	1 ρ s	∂σ s ∂x	= 0	(B.3)

s (x, t) being the velocity field, and σ s (x, t) the stress, (B.1) may be rewritten as which is equivalent with (B.1). Now let's write (B.6) in matrix form

Remerciements

of the wall. Additionally, one can also find a second peak which appears at t = 0.6 s, corresponding to the sub-figure (k) in Figure 6.20. It is the moment that the backflow impacts on the right side of the elastic wall. In Figure 6.22 one observes that the fluid kinetic energy stopped to increase at the moment t ≃ 0.2 s because of blocking effect of the elastic wall. And then, as time goes on, the fluid energy is dissipated numerically due to the use of Runge-Kutta time integrator combined with approximate Riemann solver for each pair of fluid particles.

Since no physical dissipation model is introduced in this coupled system, the energy dissipation may be only induced from three aspects: (1) numerical scheme for fluid subdomain; (2) numerical scheme for solid sub-domain; (3) the used coupling procedure. For the first aspect, as we discussed previously, the fluid sub-domain is simulated by an energy-dissipating scheme. Whereas, in Finite Element method, the explicit Newmark time integrator (β = 0, γ = 0.5) for the structure is well known as an energy-conserving numerical scheme. To know if the coupling method dissipate or inject energy at the interface, we calculate the interface energy (5.7) by (6.8).

Appendix A Analytical solution of the 1-D piston problem

The 1-D problem of mass-spring system coupled with a column of acoustic fluid is rebuilt in such a coordinates system An acoustic fluid in the 1-D tube can be governed by the wave equation

where u f (x, t) denotes the fluid displacement, c ref

where p b denotes the background pressure, ρ ref

where A s is the section area of the tube. And the initial conditions may be given as

We will apply the method of separation of variables to resolve this 2 nd order partial differential equations (A.1), with two boundary conditions (A.3) and two initial conditions (A. [START_REF] Belytschko | Element-free Galerkin methods[END_REF].

Assume that the trial solution u f (x, t) can be rewritten as the product of two functions u f (x, t) = X(x)T (t) (A.5) Substituting (A.5) in the partial differential equation (A.1) gives

where T ′′ (t) denotes the second-order derivative of the time function T (t), and X ′′ (x) the second-order derivative of the space function X(x).

If X(x) = 0 and T (t) = 0, (A.6) can be written as

where λ is a constant which is independent of x or t.

Hence we have

Let's now consider the differential equation of X(x)

which possesses a general solution

where C 1 and C 2 are two coefficients to be determined.

With the boundary condition at x = 0, we have

hence the general solution becomes

With the first initial condition in (A.31), we obtain

Through the superposition theory, we have