Guillaume

Saira

Wen Wiebke

Gilles-Alexis Martin

Hélène Magdalena Nicolas

Pascal Michael

Keywords: COUPLAGE FLUIDE-STRUCTURE, ELEMENTS FINIS, SMOOTHED PARTICLE HYDRODYNAMICS, IMPACTS, DECOMPOSITION DES DOMAINES FLUID-STRUCTURE INTERACTION, FINITE-ELEMENT ANA-LYSIS, SMOOTHED PARTICLE HYDRODYNAMICS, IMPACTS, DOMAIN DECOM- . . . . . . . . . . . 34

Je voudrais remercier tout d'abord mes encadrants, Alain Combescure et Jean-Christophe Marongiu, pour m'avoir fait partie de cette aventure où la plus grande récompense que j'en tire c'est la joie d'avoir pu travailler avec vous et surtout la grande quantité des connaissances que vous m'avez transmises pendant ces trois dernières années. Un grand merci aussi à Etienne Parkinson pour m'avoir donné la chance de faire partie de cette magnifique initiative.

Egalement, je remercie Michaël Brun, qui m'apporté ennormement de l'aide et des conseils précieux pendant ces travaux.

Je remercie chaleureusement les membres

2.4

Reference problem for establishing the solid sub-domain's governing equations .

2.5

Illustration of the outward pointing normals, n s and n f , for each sub-domain 3.1 Overview of the coupling procedure when the solid uses a smaller timestep than the fluid, i.e. ∆t f = m∆t s .

3.2

Flowchart broadly describing how the algorithm is implemented for the case when smaller time-steps are used in the solid sub-domain . . . . . .

3.3

Overview of the coupling procedure when the fluid uses a smaller timestep than the fluid, i.e. ∆t s = m∆t f and information from the interface comes only at the beginning of each macro-time-step t 0 = t n . . . . . . . 3.4 Flowchart broadly describing how the algorithm is implemented for the case when using smaller time-steps in the fluid sub-domain . . . . . . . .

Introduction

Fluid-structure interaction phenomena occurring over extremely short time scales often require expensive experiments, motivating the ongoing development of numerical tools Introduction of rock projectile impacts on turbine blades as well as other research topics at which the present project aims to shed light upon.

Hydroelectric power generation

Hydroelectricity remains one of the main sources for energy production of renewable nature. As of today, among all renewable energy sources, energy produced through hydroelectric stations represents 70% of energy generation and almost 17 % of the overall energy produced globally [REN 16]. For industrialised countries, the exploitable sources for producing hydroelectricity have reached their peak, but among developing countries, the hydroelectric presence in the energy production market is continuously rising. The total energy output coming from hydroelectric sources in 2015 was estimated at 3 940 TW h [REN 16].

The main advantages making hydroelectric energy generation an attractive and viable option are:

• Its flexibility, which stems out from the fact that a power station can be relatively simply halted and restarted in a matter of minutes following output power demand.

• The relative low cost of the energy produced [REN 16]. Hydroelectric plants can be relatively inexpensive to build, they require few personnel to operate and their expected functional lifetime can span across multiple decades.

• The possibility to stock the energy producing resource for use during peak demand periods. Additionally, water reservoirs put in place by the construction of hydroelectric stations can provide other services to reinvigorate the economy of the region such as aquaculture or by attracting visitors to the region.

Despite many of the interesting attributes that hydroelectric power generation provides, some drawbacks exist that should also be considered, such as:

• A disruption to the initial habitat in the region where a hydroelectric power station is built. Usually through the flooding of the regions to build a dam, many of the flora and fauna as well as the local human inhabitants must be displaced. This can have enduring consequences on the biological equilibrium of the region.

• The accumulation of silt within the premises of the station and reservoir. Silt particles transported downriver by the flow of rivers and streams can cause the station to become non-operational if the reservoir eventually becomes saturated.

As with any source for energy generation, hydroelectricity production presents both advantages and drawbacks. The goal of the work in the present study aims at providing a reliable method for numerically simulating potential problems that might arise throughout the exploitation of a hydroelectric power station thus contributing in the long run to providing safer conditions for hydroelectric energy generation.

The Pelton wheel

The project behind the current works aims at characterizing erosion damage undergone by hydraulic turbines of the Pelton type. These type of turbines have been around since their inception by the American inventor and entrepreneur Lester Allan Pelton in the late 19th century [WIL 74]. Opposed to the centuries old water wheel, where most of the energy was transferred mainly by the dead weight of water on each of the partially submerged blades or buckets, the design proposed by Pelton allowed to transfer the momentum energy from the water impinged at higher velocities on each of the buckets through a converging nozzle. The design of the buckets has also evolved over time in order for water to transfer as much momentum as possible to the turbine. In fact water in contact with the bucket, performs a U-turn and exits the bucket at a much lower velocity, thus Pelton turbine-based generators are single stage.

Compared to other hydroelectric turbines commonly used today, such as turbines of the Francis or Kaplan type, Pelton turbines are mainly employed in contexts where a high dynamic heads (greater than 100 m) exist with lower flow rates (discharges lower than 50 m 3 /s). Pelton turbines come in a variety of sizes ranging from, large ones capable of outputting 440 MW to small ones that are just a few centimeters wide and handle low flows.

In terms of design specifications, the magnitude of the tangential velocity of the turbine is usually set to be half of the jet's velocity impinging the bucket [WIL 74] Introduction With g being the acceleration due to the gravity field and h is the dynamic head. For a wheel with a diameter D w , the angular velocity of the turbine is set at:

ω = V t D w /2
For the current study, a high dynamic head will be considered, h = 2000 m. Hence, the velocity of impinging jet will be V j = 200 m/s and the tangential velocity of the wheel will be V t = 100 m/s. Thus, the relative velocity at which the jet will be impacting the turbine blade or bucket is V r j = 100 m/s assuming the jet and tangential velocity vectors to be collinear.

The PrEDHyMa project

The PrEDHyMa project, whose acronym stands for Predicion of Erosion Damage on Hydraulic Machines, aims at studying four possible sources of damage that hydroelectric turbines undergo during their functional lifetimes. The project is aimed particularly at characterizing damage on Pelton turbines which are used to produce electricity when high heads are available notably in mountainous regions. The four sources on which the project focuses are:

• Gravels and stones impacts: The impact of gravels of stones found in streams during the turbine's operation can lead to strong damage in its buckets. The possibility of this event happening is not very high since precautions are normally taken in order to prevent its occurrence, nevertheless, the simulation of this event is important since the violence of the impact can heavily damage the turbine.

• Droplet impacts: The impact of water droplets on the surface of the blades of a hydroelectric turbine can lead to damage through repeated exposure to this type of phenomenon. The damaging mechanism is not the same as in the previous case since erosion on the surface doesn't arrive suddenly but through prolonged exposure.

• Hydro-abrasive erosion: The presence of small sediments in streams and rivers provides a source of erosion that is often responsible for the damaging of Pelton turbines. Due to the small size of the solid sediments, visible damage usually arrives after prolonged exposure to this type of event. The phenomenon studied within the scope of the project take place over different time scales. The first phenomenon, due to the size of the impacting body can cause immediate visible damage on the bucket of a turbine whereas damage induced by the other three phenomena might become apparent through repeated occurrences thus a fatigue failure prediction tool is required to evaluate its extent. The current work is focusing on the former type of phenomenon so only a fluid-solid solver coupling tool capable of simulating impacts between two structures in a multi-time step configuration will be developed. A tool capable of using the fluid-solid coupling capabilities of the current work to estimate damage by repeated cycles will also be developed within this project but not used in the current case [COU 15].

Fluid and solid solvers

The current thesis will focus entirely on building a numerical tool for coupling specialized fluid and solid solvers in order to simulate the impact of gravels and stones on a Pelton turbine bucket turning at full regime. These simulations will allow to evaluate the extent of the damage sustained by the target structure if such event were ever to occur.

The specialized solver for the solid sub-domain is called Europlexus is developed by Introduction the French CEA (Commisariat à l'Energie Atomique ) and the European Commission's Joint Research Centre (JRC). This solver, mainly used to carry out non-linear transient dynamic calculations, is mostly based on the Finite-Element Method and the use of a second order time integration scheme [START_REF]Commisariat á l'énergie atomique[END_REF]. Example applications of this solver include [FAU 03, CAS 03, PRO 16].

The specialized fluid solver that will be used for this study is called Asphodel and is developed by Andritz Hydro. This code uses an explicit time integration scheme and is based on the mesh-less technique known as the Smoothed Particle Hydrodynamics (SPH) and allows to simulate free surface flows interacting with complex solid geometries without the need of using sophisticated interface determination techniques. Details on the development of this solver as well as example applications can be found in [ MAR 07,NEU 14,REN 15].

The coupling of the solvers described here will be the at the core of the research proposed by the present work. As it will be explained in detail in the current text, the use of explicit time-integrators for each sub-domain constraints greatly the time-step size of each solver an thus the total duration of the simulations may be quite large. In order to give a certain degree of independence to the choice of the time-step on each solver, a multitime step algorithm is proposed. The coupling of the solvers involved has been done using different communication solutions between the solvers and a coupling software (FIFO pipes, binary files). While interesting and somewhat complex, the details regarding these technical aspects will be omitted in order to avoid misguiding the reader from the real research subject proposed by the current work. More details regarding the communication tools used to couple the solvers in this thesis can be found in [FOL 11, WIL 05, SNI 98].

Development of a multi-time step fluid-structure coupling strategy

The current work focuses on building upon a coupling strategy to simulate transient fluidstructure interaction. A first approach of the method was developed by [ ZHE 13], [LI 15]: it couples a solid Finite-Element and a fluid SPH (Smoothed Particle Hydrodynamics) solvers in a non-intrusive and synchronized manner. By imposing the normal velocity continuity at the interface, this coupling method ensured that neither the energy production nor dissipation occurred at the interface, thus guaranteeing the coupling simulation's stability over time. The codes use different time integrators for each sub-domain: 2nd order Runge-Kutta scheme for the fluid and an explicit Newmark time integrator for the solid. The method has been developed for applications using the same time-step in all sub-Introduction domains. This work presents a method capable of handling different time-steps in each sub-domain. It is based on a dual Schur method, called GC method [GRA 01], which extends the widely used FETI techniques [FAR 91]. The GC method, initially proposed for structural dynamics, is able to couple any time integrators from the Newmark family with their own time-step depending on the sub-domain under consideration. The method is extended here to the coupling of a Newmark time integrator and a Runge Kutta scheme with large interface displacements.

There are basically two ways regarding how FETI-based techniques can be used to couple sub-domains with different time scales. In the first group of algorithms, the kinematic condition at the interface is satisfied at the largest time-scale, such as in [PRA 04] and [MAH 11]. The second type establishes that kinematic continuity should be prescribed at the smallest of time scales, presented by [GRA 01]. More recently, Brun et al. [BRU 15] proposed a general framework for building energy conserving methods at the two time scales (kinematic continuity at the large or the fine time-steps), by cancelling the interface pseudo-energy involved in the so-called energy method, this was originally proposed by Hughes for proving the stability of hybrid explicit-implicit time integrators [HUG 12]. In all these methods it has been proved through the energy method, that kinematic continuity must be prescribed in terms of velocity to guarantee the stability of the coupling algorithms.

The aforementioned coupling methods are endowed with advantages and drawbacks as, for example, the macro based methods turned out to be more suited for parallel computing. However, since the interface constraint condition is only verified at the macro time-scale, an accurate tracking of the interface motion becomes more challenging. The micro-based multi-time-step coupling techniques handle more easily strong geometrical and material non-linearities at the fine time scale. However, some small amount of numerical dissipation can occur at the interface, due to the required interpolations of quantities at the large time scale. The present work targets the simulation of highly non-linear phenomena, both from the geometrical and the material perspectives, occurring in fluidstructure interaction problems. The GC method, due to its micro time-scale coupling characteristics, has thus been chosen.

As previously explained, the sub-domain computations are performed with external specialized solvers. These solvers are linked by an external coupling software, handling only the interface computation derived from the normal velocity continuity at the interface between the fluid and solid sub-domains. The unknowns of the interface problem are the Lagrange multipliers, which represent the interface pressure. They are computed Introduction by the external coupling software and sent to the sub-domains so as to complete their computations for a given time-step.

This thesis work is organized as follows: The first chapter provides a state of the art regarding techniques to couple solvers, the different approaches to modelling fluidstructure interaction, domain decomposition methods and solid contact modelling. The second chapter presents the governing equations for each sub-domain and the continuity equations at the interface as well as the discrete numerical schemes used to numerically describe these equations in time and space. The third chapter introduces the multi timestep coupling strategy, which is the main contribution proposed by this work. Finally numerical examples allowing to validate the proposed technique and state its limits are given in the last chapter along with impact simulations that give an insight to the damage undergone by turbine buckets when hit by rock projectiles carried by fluid flows.

Chapter 1

Domain decomposition methods for fluid-structure interaction problems

The purpose of this chapter is to give the reader an overview of the numerical tools that can be applied for fluid-structure interaction problems allowing for temporal and spatial decomposition of the numerical domain

Introduction

The simulation of the phenomenon presented in the introduction requires the use of techniques allowing to numerically couple heterogeneous physical domains. Domain decomposition methods exist that allow to carry out simulations with different spatial and temporal decompositions and have been the subject of ongoing research since the dawn of the development of numerical analysis tools [SCH 69, BEL 78].

In the present day, we can find a large number of techniques allowing to effectively solve a large and complex problem as the superposition of smaller and simpler ones [FAR 91, GRA 01, PRA 04]. Such techniques not only may take advantage of today's massively parallel computer architectures but allow to model a seemingly complex numerical problem as a collection of simpler ones while more readily adapting to the specificities of each simulation.

In the current work we will focus on giving the reader an insight into the main types of coupling tools that allowed us to numerically couple different physical domains using heterogeneous time scales. To do so we will separately present:

(ii) Dynamic equilibrium across the interface:

F 1 Γ + F 2 Γ = 0
By meeting the above statements as best as numerically possible, one can produce stable and accurate coupling algorithms allowing to link different physical media while taking advantage of the calculators parallel architecture through efficient domain decomposition.

Domain decomposition methods for fluid-structure interaction problems

Coupling of fluid-structure solvers

Numerical tools used to simulate the interaction of different physical fields allow -in many cases-to obtain more detailed and precise solutions to problems with respect to simulations in which hypotheses and model reductions can lead to neglecting the interaction of a physical system with surrounding physical fields. However, one should take into account the very logical fact that by adding the influence of other physical fields to the initial problem, the problem becomes increasingly more complex to solve as not only one must ensure that the additional equations induced by the existence of new fields to be satisfied, but a series of additional conditions must be met at the interface separating the physical systems so as to ensure the accuracy of the solution and stability of the new problem . This latter aspect is of utmost importance and can become a determinant factor when deciding if a physical system should be modelled by taking into account the influence of surrounding physical fields or not. If the former were the case, two major ways of dealing with how the dynamic coupled problem can be envisaged [FEL 01] 1. Monolithic approach : The whole physical system is treated by a single solver including the problem at the interface between sub-domains and all of its entities are advanced in time synchronously (1.1 (a)).

2. Partitioned approach : Each physical entity is handled by a separate solver or solvers, the problem at the interface may need to be solved in order for time to be advanced separately (1.1 (b)).

As Figure 1.1 (a) broadly illustrates, a monolithic treatment to solve a numerical problem, involves generally a single solver -which can be running over a single or multiple calculation cores -that handles all the sub-domains involved in the coupled calculation.

With this type of approach the equations governing each physical sub-domain (Ω 1 and Ω 2 ) and the conditions at their shared interface (Γ) are discretized and the new configuration is obtained in a synchronous manner. In contrast, 1.1 (b) depicts the general manner in which a partitioned approach handles the access to the updated configuration for the physical system Ω(t + ∆t). Using this type of treatment, each sub-domain present is updated in time independently and continuity conditions at the interface may be treated separately. In order for a partitioned approach to remain conservative, special care must be taken when solving the coupled equations so as to ensure that the continuity conditions are respected as best as numerically possible [FAR 91, GRA 01, PRA 04, LI 15].

Γ Ω

1 Ω 2 Ω 2 Γ Ω 1 Ω(t) Ω(t+Δt) (a) Ω 1 Ω 2 Ω 1 Ω 2 Γ Γ Ω(t) Ω(t+Δt) (b)
Figure 1.1: Illustration of a numerical problem in which physical domain Ω is decomposed into sub-domains Ω 1 and Ω 2 with boundary Γ.Treatment of the numerical problem through a monolithic approach 1.1 (a) and a partitioned approach 1.1 (b)

In general, when solving a problem through a partitioned approach, a prediction is made on one of the fields' variables which allows to find the solution of the whole physical system in a staggered manner (Figure 1.2). The accuracy of the solution depends naturally on the quality of the predicted quantities. On the other hand, when solving a problem through a monolithic approach, the whole system of equations is solved directly (through direct or iterative algorithms), allowing for all the fields' variables to be updated at once.

Figure 1.2 shows a broad example in which a physical system is updated staggeringly using a partitioned approach with synchronization at the interface. In this example, the new configuration of sub-domain Ω 1 is estimated and used to obtain (1) an estimated solution of the interface continuity equations. Then the predicted interface status is sent to sub-domain Ω 2 (2) that uses it to find its new configuration (3). Next, Ω 2 new configuration is used (4) to verify the continuity equations giving access to the updated configuration of the interface. A convergence test can be done at this state to make sure that the continuity equations are satisfied at the interface which can lead to iterations between Ω 2 and Γ. After satisfying the interface continuity equations, the updated interface status is sent to sub-domain Ω 1 (5). Using this information, Ω 1 can be updated (6) and the process can be repeated.

The advantage of one type of coupling treatment over another are dependent on the nature of the problem as pointed out in [FEL 01, BEL 78]. One may be inclined to treat the simulation of a physical system in a monolithic way if the simulated fields share comparable characteristic time lengths and time scales or if the interaction effects between the physical fields are widespread within the domain. On the other hand, it can be advisable to

Ω 1 Ω 1 Ω 2 Ω 2 Γ 1 2 3 4 5 6
Figure 1.2: Partitioned approach typical sequential algorithm with synchronization treatment at the interface Γ between sub-domains Ω 1 and Ω 2 try a partitioned approach if the problem features localized interaction between physical fields or the time-scales characterizing the involved fields are quite different.

Regarding implementation convenience, each kind of treatment has both its share of advantages and inconveniences, as the monolithic approach allows to shorten development times and can facilitate software maintenance, however, computational efficiency can be greatly hampered as a result of the large linear systems the single solver might need to handle by itself.

Partitioned approaches are desirable when external specialized solvers are available and only the way in which they interact among them is to be added. This allows for a higher degree of independence/flexibility for each specialized solver and foments software modularity (e.g. addition of turbulence models for a fluid solver or erosion models for a solid solver if considering fluid-structure interaction). Nevertheless, in order to maintain acceptable levels of accuracy and stability, the solution algorithm needs to be well thought-out and implemented. This last statement branches out from the fact that when solving a coupled system with a partitioned algorithm, a time lag might exist between the actual moment in which each sub-domain is updated to a new configuration. If precautions are not taken to minimize the effect of this configuration delay between physical fields, the algorithm might be non-conservative and might diverge or give spurious results. For this reason, it is important to translate as best as possible the kinematic continuity and dynamic equilibrium conditions brought-up in the introduction to the current section.

In our case, we are considering the use of specialized solvers for each sub-domain and an external coupling software to be the link between them. From this description and as Figure 1.1) conveys, we can advance that a partitioned approach seems as the most convenient way to carry out our simulations. Thus a synchronized partitioned approach Fluid-structure interaction -the different approaches is to be used allowing to make use of different time scales in each sub-domain. This approach will be explained in more detail in the following sections, the essential thing being that by synchronizing the way in which the kinematic quantities are obtained for each solver, a stable and accurate solution algorithm can be established even if the time integration is done separately for each sub-domain.

Fluid-structure interaction -the different approaches

After a brief introduction on the coupling approaches that can be used for simulating physical systems, we will now focus on the core subject of the current thesis by taking a look at some methods that exist for simulating physical systems where the effect of fluidstructure interaction can be considered preponderant. Due to vastness of the subject, this section cannot be exhaustive, but will focus on giving the main details regarding some of the most well-known techniques currently in use and will attempt at justifying the choices made for the simulations that were carried out within the scope of this work.

The approaches put forward in the current section can be regrouped into two larger groups : Analytical methods and semi-discrete methods. The common point regarding all of these methods is that they can all be used to solve FSI problems while respecting the continuity conditions evoked earlier, however, as we will see, their suitability can vary with respect to the different factors characterizing the FSI problem at hand such as the size, the geometries and the nature of the fluid flow and the solid's response.

Analytical Methods for fluid-structure interaction

Analytical methods have long been employed for solving fluid-structure interaction problems. In fact for cases where the solid and fluid geometries are quite simple and several assumptions can be made about the flows and the amplitude of the structure's response, these methods can be second to none when considering their ease of implementation and reliability. The works of [PAI 98] and [DEL 01, BEL 80] highlight the main aspects behind these methods and the cases in which their use is advisable are put forward.

Boundary Integral Methods

Boundary integral methods allow to find the forces and the moments applied to completely or partially immersed solid volumes within potential flows (inviscid and incompressible 1. Domain decomposition methods for fluid-structure interaction problems fluid) or acoustic flows (compressible fluid). By using simplification tools such as the divergence theorem or the Reynolds transport theorem, the resolution of the partial differential equations gives access to the induced forces on the boundary between the fluid and solid. Furthermore, Green's identities [STR 92] can be used and the FSI problems can be further simplified for cases where Green functions exist which include potential or acoustic flows (1.3). The greatest advantage of boundary integral methods, as the name suggests is that the calculations are done on the boundary only, thus reducing the dimension of the problem by one. Alternatively, for incompressible and acoustic flows, for which Green's functions exist as solutions to wave equations, using the divergence theorem to obtain Green's identities, the problem is once again solved only on the interface and not on the bulk. If two harmonic scalar functions, φ and ψ, are considered, Green's second identity can be written as:

S φ ∇ψ • n dS = S ψ ∇φ • n dS (1.1)
where n is the normal vector pointing outward from the volume surrounded by the surface S.

Replacing ψ by the Green function which is a solution to the wave equation and φ by the pressure or velocity field, and applying the appropriate boundary conditions for the pressure and the velocity at the interface S, the value of the pressure and velocity fields can be found anywhere within the fluid sub-domain even in the presence of field sources and sinks within the fluid sub-domain [WAN 08].

Furthermore, this method can be extended to complex geometries by discretizing the whole fluid-structure interface [SIG 15]. This is the underlying idea behind the Boundary Element Method which can be quite efficient for finding the solution of fluid-structure problems with equal accuracy than the FEM method but with the advantage of being solved on a support domain with a lower overall dimension than the original problem. Nevertheless, flow configurations for which Green's functions exist as solutions to the governing differential equations are limited to cases of potential or acoustic flows, hence for any other type of flow, other analytical or numerical techniques are more appropriate.

Dimensional Analysis for Fluid-Structure Interaction

Dimensional Analysis, which is a conventional technique to study problems in fluid mechanics, can be extended to be applied in fluid-structure interaction problems as exposed in the works by [DEL 01]. The underlying idea behind this method used to study fluid-structure interaction problems is to reduce the total amount of parameters and variables involved the dimensionalized equations to the strict minimum using the Vaschy-

Buckingham-π theorem [MAC 71].
Using this method, the conservation equations in the solid and the fluid sub-domains as well as continuity equations along the interface can be expressed as a function of dimensionless numbers.

Regarding commonly used dimensional parameters in analytical fluid-structure interaction we can list:

• For both domains, x is the position vector of the current point, t is the time parameter, L is a characteristic length of the domain under consideration and g is the norm of the gravity vector

• For the fluid domain, µ is the dynamic viscosity, c 0 is the reference sound speed, U 0 is a reference velocity within the domain, p 0 is the reference pressure and ρ 0 F is the reference density of the fluid • For the solid domain, E is the material's Young modulus, µ is the Poisson's coefficient, ξ 0 is the reference displacement of the structure and ρ 0 S gL is the reference 1. Domain decomposition methods for fluid-structure interaction problems density of the solid material.

Using the cited dimensional parameters, the following dimensionless quantities can be derived for both the fluid and the solid:

• x/L is the dimensionless position • t/T f luid or t/T solid (with T f luid = L/U 0 and T solid = L/ E/ρ 0 S )
is the reduced velocity and is used to compare the time scale of the studied phenomenon with respect to the time it takes for a point in the fluid continuum to travel or propagate through the characteristic length L

• Π = ρ 0 F U 2 0
p 0 is the inverse Euler number and it is used to express the relationship between a local pressure variation and the kinetic energy per volume

• R e = ρ 0 F U 0 L µ
is the Reynolds numbers which quantifies the importance of convective flows with respcect to viscous forces

• F R = U 0 √
Lg is the Froude number which measures the importance of inertial and gravitational forces

• M = U 0
c 0 is the Mach number which measures the compressibility effect

• D = ξ 0 /L is the reduced solid displacement and measures the effect of large displacements

• G = ρ 0 S gL
E measures deformation on the structure by the simple effect of gravity. This number can be replaced by the reduced velocity U 0 /c s which compares the fluid flow velocity U 0 to the elastic wave propagation speed in the solid c s .

.

In order to couple both domains, one last dimensionless parameter, A, should be considered which takes account the coupled nature of the problem. According to the nature of the fluid-structure interaction problem considered, A can take different forms such as: From the description of the parameters , it can be seen that certain parameters describe the fluid flow (Π, R e , F R , M) while others describe the solid motion and behavior (µ, D, G) and A (M , C Y ) describes the interaction between the solid and fluid domains By writing the conservation equations along with the interface continuity equations in their dimensionless form, an analysis of the relevance or negligibility of certain dimensionless parameters can be carried out when knowing the fluid's and solid's physical characteristics as well as the flow conditions and the nature of the solid response. This dimensional analysis can help to greatly simplify equations and allowing them to be solved analytically.

Using dimensional analysis, different kinds of fluid-structure interaction problems can be solved using hypotheses based on the physical conditions governing the problem. For the case of the solid floating on a fluid at rest (Figure 1.4 (a)), the reference fluid velocity can be neglected (U 0 = 0) and the displacement of the solid (ξ = q(t)Φ(x)) can be simplified to depend only on time (rigid body motion). These hypotheses can lead to establishing a governing equation for the solid motion in which the effect of the fluid appears only as an added stiffness:

q(t) + M F 2 D q(t) = 0 (1.2)
Regarding the case illustrated in Figure 1.4 (b), the flow velocity (U 0 ) can be considered preponderant with respect to the motion or vibrations in solid. If we consider only the rotation of the airfoil, after some simplifications, the following governing equation of motion can be established: 1. Domain decomposition methods for fluid-structure interaction problems

θ(t) + (1 -2πC Y x L )θ(t) = 0 (1.3)
Hence, once again, an added stiffness is applied to the system by means of the fluidstructure interaction. For this case, however, the sign of the stiffness can be change as a function of the flow velocity and the position of the barycentre of applied fluid forces on the airfoil (x). By studying the sign of the added stiffness, configurations for which the flow-induced vibrations grow exponentially can be found thus giving an insight to flow-induced instabilities on the solid.

The interest of carrying out dimensional analysis for fluid-structure interaction problems comes mainly from the fact that a long range of problems can be treated if the appropriate hypotheses are taken in order to establish the governing equations. Nevertheless, the existence of these simplifying assumptions many times requires an a-priori knowledge of the flow conditions and the associated structure's response. For cases in which fluid flow has reached its steady-state regime, these hypotheses might be more easily obtained and the use of this method can be straightforward. However, for transient, free surface, fluid-structure interaction problems, the establishment of the governing equations can be more delicate if at all possible.

Method of characteristics

The method of characteristics is somewhat similar to the method of the retarded integrals for solving the boundary integral problem. For this method, the solution to hyperbolic PDE's is found along the characteristics of the wave-type solutions, allowing, on one hand, to obtain a simpler system of equations to solve based on a ODE (since the problem is only solved through the characteristics, a intrinsic relationship between the spatial and temporal variables can be established), and on the other hand, to extrapolate the value of certain fields from within the fluid sub-domain into the solid interface and vice-versa. The method of characteristics is mostly applied of the simulation of compressible flows where the presence of shocks is preponderant. However, its implementation can be a obstacle when dealing with complex geometries and free-surface flows (as opposed to using discretization techniques such as FEM or FVM). The reader can find more information regarding this technique by consulting [JOH 82].

Semi-discretization methods for fluid-structure interaction

Semi-discretization methods offer the greatest versatility regarding the implementation of solution algorithms for FSI problems. These methods allow, by discretizing in space the governing PDE's on each sub-domain, to obtain an approximate ODE, which can be integrated in time. Semidscretization methods are the one of the most encountered tools used for solving FSI problems nowadays. These methods include the Finite-Element method, the Finite-Difference method, etc.. and can be coupled in such a way that the optimal method for each sub-domain can be used independently of the semidiscretization methods used for the neighbouring sub-domains. However, these methods might require, in certain cases, the use of finely discritized zones, in order to increase the accuracy of the solution. This causes a surge in the number of degrees of freedom (DOF) of the problem, leading to a higher demand of computational resources in order to obtain the solution to the coupled problem, particularly in the presence of implicit solution algorithms. Over the past years, numerical techniques have been developed in order to allow obtain more precise solutions, even when using a coarsely discretized meshes [ABG 03]. These and the ongoing progress in computer science, allow for these type of methods to be among the most popular choices when implementing an FSI solution algorithm. For these reasons, running simulations using semidiscretized methods seems as the most logical choice regarding the solving of the numerical problem aimed within the scope of this thesis. In the next section, we will focus on giving some details regarding some of the families of semidiscretized methods used today and a justification will be given regarding the choice of the semi-discretized method used to carry out our simulations.

The main objective within the scope of this work is to carry out simulations using semi-discretized methods for fluid-structure interaction. In this section, the main methods used to solve the fluid and solid sub-domain equations will be introduced in a broad manner. The mathematical formalism expliciting how to obtain the semi-discrete system of equations from the governing equations for each sub-domain will not be explained for each method introduced. This will only be done for the chosen methods in the following chapter of this thesis.

Finite Difference Method

One of the first widely-spread discretization techniques used to analyse the behaviour of structures and fluids among with their interaction was the Finite Difference Method (FDM). This technique's success relies on the simplicity of passing from a continuous 1. Domain decomposition methods for fluid-structure interaction problems PDE in its strong form to a system of discretized ODE's through a Taylor expansion. By using numerical enhancements such as multi-grid techniques [VEN 00, MOR 05], the solution procedure can be sped-up even when dealing with many degrees of freedom

Ω Ω d Γ u| Γ = u 0 u(t=0) = u t0 i, j i-1, j i+1, j i, j-1 i, j+1 h i h j x y x(i) y(j)

PDE

Figure 1.5: Schematic illustration of the way in which a continous domain is approximated using finite differences This type of method is simple to derive and to implement for the discretization of the fluid and the solid conservation equations as well as for the continuity conditions along the interface. Nevertheless, this method is aimed to discretize physical domains which would follow ordered mesh topologies since the imposing of boundary conditions can become difficult when dealing with geometrically complex domains. In order to use this method with free-surface flows, this might require the use of complex numerical tools in order to efficiently track the fluid-solid interface or the fluid's free-surface topological definition. For these reasons the Finite Difference method cannot be used to carry out the free surface flow simulations that the current work necessitates. More information on this method can be found in [GRO 07].

Finite Element Method

To overcome some of the drawbacks of the Finite Difference method, the Finite Element Method (FEM) was introduced later on [HRE 41, ZIE 77]. A finite element method discretization is based upon a piecewise representation of the solution in terms of specified basis functions. The computational domain is divided up into smaller domains (finite elements) and the solution in each element is constructed from the basis functions. The actual equations that are solved are typically obtained by restating the conservation equation in weak form: the field variables are written in terms of the basis functions, the equation is multiplied by appropriate test functions, and then integrated over an element. Since the FEM solution is in terms of specific basis functions, a great deal more is known about the solution than for either FDM. Thus, the choice of basis functions is very important and boundary conditions may be more difficult to formulate. Again, a system of equations is obtained (usually for nodal values) that must be solved to obtain a solution. The Finite Element Method offers the possibility of numerically solving problems with complicated meshes and boundary conditions. For this reason it becomes the numerical tool of choice of a great number of researchers seeking to model the behaviour of structures even when undergoing large displacements. For fluid mechanics problems, the use of a Lagrangian mesh would be desirable in order to be able to track free surfaces or fluid-structure interfaces. However this would prove to be quite a difficult task considering the large distortions such a mesh would undergo for certain types of fluid flows. For the former reason we have chosen to carry out the numerical simulations within this thesis work making use of the Finite Element for the solid sub-domain. Concerning the fluid sub-domain, two numerical methods which can be considered better suited to be used for numerical simulations will be briefly presented next.

Finite Volume Method

A finite volume method (FVM) [LEV 02] discretization is based upon an integral form of the PDE to be solved (e.g. conservation of mass, momentum, or energy). The PDE is written in a form which can be solved for a given finite volume (or cell). The computa-1. Domain decomposition methods for fluid-structure interaction problems tional domain is discretized into finite volumes and then for every volume the governing equations are solved. The resulting system of equations usually involves fluxes of the conserved variable, and thus the calculation of fluxes is very important in FVM. The basic advantage of this method over FDM is it does not require the use of structured grids, and the effort to convert the given mesh into structured numerical grid internally is completely avoided. As with FDM, the resulting approximate solution is a discrete, but the variables are typically placed at cell centers rather than at nodal points. This is not always true, as there are also vertex-centered finite volume methods. In any case, the values of field variables at non-storage locations (e.g. vertices) are obtained using interpolation. Due to its simplicity (implementation is similar to the FDM) and its adaptability when it comes to the simulation of complex fluid flows with unstructured geometries, the Finite Volume Method is one of the most widely used techniques for the numerical simulation of fluids. Despite these advantages, the Finite Volume Method presents the same shortcomings as all other mesh-based method which is mainly its inability to handle large mesh deformations or the difficulty it represents to effectively track interfaces between the fluid and other physical media. Hence, a meshless method will be introduced next, as we believe represents the best alternative for simulating the fluid flows considered within the scope of this work.

Smoothed Particle Hydrodynamics method

The Smoothed Particle Hydrodynamics (SPH) [GIN 77, LUC 77, VIL 99] is a numerical method used to solve partial differential equations on a cloud of calculation points x i in a computational domain Ω ⊂ R d . It is a meshfree method because there is no connectivity between the calculation points. Traditionally, the SPH operators, that will be introduced below, are used to discretize the Euler or Navier-Stokes equations in Lagrangian form.

Fluid-structure interaction -the different approaches

The calculation points are interpreted as particles that follow the flow in Lagrangian motion. Additionally this method holds several similarities with the Finite Volume method presented previously. In fact, for Euler or Navier equations both methods need to compute the time rate of change of the flow variables averaged over a control volume by summing numerical fluxes over neighbouring cells or particles. However for the SPH method, neighbours to a particle i change constantly and are located within a finite distance depending on the size of a weight or kernel function. Thus the numerical stencil of the SPH method is much larger than that of the Finite Volume method which only takes into account its adjacent neighbours.

i The SPH method has been found to be specially well suited for capturing dynamic free surface flows. The fluid particles follow the flow and they are able to capture the position of the free surface. In addition, no mesh has to be created to connect the calculation points that are displaced with the flow. Additionally, it is possible to use a variant of the SPH method, named the SPH-ALE Method [VIL 99, MAR 07], which allows further flexibility for imposing boundary conditions. For these reasons, the SPH method seems to be the best suited choice for the carrying out the calculations targeted within this work. On the one side, free surface flows are naturally taken into account by the method and on the other hand, the method can be made compatible for interaction with complex solid geometries.

One last interesting feature regarding the simulation a fluid driven stone impacting a solid target is that SPH particles can adapt easily to the relative motion of the projectile within the flow while using a mesh-based technique this will be computationally expensive if at all possible.

1. Domain decomposition methods for fluid-structure interaction problems 1. 3

.2.5 Conclusion on the semi-discretization methods

Over the course of this subsection, a very succinct overview of some of the semidiscretized methods that can be used for fluid-structure interaction has been given. For each method, the most salient characteristics have been reviewed and a discussion about the suitability of each method to handle the simulations that this work demands has been done. As already pointed out, the Finite-Element method seems to be the best choice to carry out simulations of the solid sub-domain due mainly to its robustness and compatibility with the modelling of complex geometries undergoing large displacements. The SPH method could also have been a possibility for the semi-discretization of the solid domain since it offers interesting characteristics allowing to more readily simulate crack propagation and solid fragmentation. Nevertheless, the finite-element method is far more mature and provides a more stable and reliable numerical framework for the simulation of solids and solid impacts with respect to what the SPH method currently offers [SWE 95, BEL 98]. For these reasons it has been decided to exclusively use the former method for the numerical representation of the solid structures in the current thesis.

On the fluid side, however, it has been discussed that the SPH method is the only numerical method that can provide simultaneously with simple free surface tracking and the ability to interact with complex solid geometries and motions. Hence, this will be the method of choice for fluid simulation in the current work.

In the following sections, more details will be given regarding the SPH and the Finite Element method and the mathematical and numerical techniques allowing to make coupled Fluid-Structure interaction calculations with these two numerical methods regardless of their differences in spatial and temporal integration schemes.

1. 3.3 Coupling of semi-discretized methods for fluid-structure interaction FSI numerical coupling techniques exist to link the aforementioned numerical methods for fluids and solids. For mesh-based numerical methods (FDM, FEM, FVM), when using Lagrangian meshes for the simulation of both sub-domains, the distortion of the mesh on the fluid sub-domain can lead to errors in accuracy as well as stability issues coming from an extreme distortion of the fluid mesh in order to follow the fluid-solid interface or to capture the motion of the free surface in a sloshing problem. This may lead code developers and users to want to impose a fixed mesh (Eulerian mesh) for the discretization of the fluid sub-domain complemented with interface-tracking techniques to follow the fluid-solid interface or the free surface position. Another solution might imply using an Arbitrary-Lagrangian-Eulerian (ALE) mesh for the fluid sub-domain. Using this type of mesh, fluid mesh nodes along the interfaces can move following the motion of the solid interface nodes which don't necessarily follow the exact motion of the fluid material points. Hence, the motion of the fluid mesh is considered arbitrary with respect to the motion of the fluid material points within the considered fluid domain. In case of large solid motions, recursive remeshing of the fluid structure interface can be necessary in order to provide stable calculation. In this brief sub-section, coupling techniques will be presented based on the description that is adopted to describe the motion of the fluid mesh within the coupled system (A Lagrangian mesh description will be considered as the default choice for the solid). This classification will allow to better understand the attributes offered by each type of mesh description and will highlight the shortcomings of using a meshbased description for the fluid in the case of the FSI simulations carried out within the scope of this work.

Eulerian mesh description for the fluid sub-domain

Eulerian mesh description for the fluid sub-domain can offer an efficient solution to simulating phenomena involving fluid-structure interaction with negligible fluid-solid interface displacement [COM 80]. As it is widely known, when adopting an Eulerian description, the mesh elements are fixed in space while the continuum can flow freely through the mesh.

In order to carry out fully coupled fluid-structure interaction simulations presenting free surface flows, an Eulerian description of the fluid mesh can be employed but interface tracking techniques need to be employed as well. One such interface tracking technique is the one presented in [PES 02] and it is known as the Immersed Boundary Method (IBM). Using this method, the fluid is discretized in space using an Eulerian Finite Difference mesh. The solid is completely immersed and it is represented by a non-conforming (with respect to the Eulerian fluid mesh) set of boundary points that are interconnected (only the boundary is represented). The solid boundary points, which follow a Lagrangian description, apply local body forces to the fluid at their respective positions and kinematic constraints involving the solid velocity and the interpolated value of the fluid velocity at the same points in space allowing to establish a strong link between Eulerian and Lagrangian variables giving access to the calculation of fields within the fluid sub-domain and 1. Domain decomposition methods for fluid-structure interaction problems the updated configuration of the solid sub-domain. In a similar manner to the Immersed Boundary method, the Ficticious Domain method [GLO 97] allows to replace the set of points applying local forces by a full volume that can be a rigid solid or can represent a solid continuum discretized using the Finite Element Method. In addition, the coupling is done by using Lagrange multipliers in order to ensure dynamic and kinematic continuity at the interface between the fluid and the solid sub-domains thus offering a way to strongly couple both physical domains. At a first glance, carrying out fluid-structure interaction simulations while using an Eulerian mesh for the fluid can seem contradictory considering that the fluid-solid interface can constantly change during a transient calculation of this nature. However, from the preceding descriptions, it can be understood that it is actually possible to carry out this type of simulations if appropriate interface tracking techniques are implemented. These techniques can either rely on the overlapping of a solid Lagrangian mesh on the Eulerian fluid mesh and imposing interface continuity with Lagrange multipliers or using volume fractions or topological functions in order to estimate the position of the interfaces in order to update the boundary conditions on the fluid mesh at each iteration. These type of techniques can be of great interest for problems in which the structures are fully immersed within the fluid and for which the movement of the solid within its fluid entourage does not involve large displacement and rotations.

ALE mesh description for the fluid sub-domain

As stated in the introduction to this subsection, an ALE fluid mesh allows nodes or cells within the mesh to move independently with respect to the fluid continuum. Using an ALE mesh description for the fluid offers significant versatility regarding fluid-structure interaction applications. As an example one can consider a NACA-type airfoil within an air flow with a given upstream flow velocity (Fig. 1.9). If the airfoil were to remain stationary during the whole simulation, an Eulerian fluid mesh would be sufficient to carry out a simulation of the flow around the structure. However if the airfoil were to change its orientation by revolving around an axis of rotation midway through the simulation, an ALE description of the fluid sub-domain's mesh would allow the mesh points and cells located near the moving interface to follow the solid wall's motion which has no a-priori reason to follow the flow of the fluid. Hence, the fluid mesh can accurately capture the movement of the structure without requiring to practice the interface tracking techniques discussed previously. Nevertheless, one major drawback concerning this description can happen if the motion of the mesh following the free surface or the fluid-solid interface induces a large distortion on the fluid mesh. If such is the case, there is necessity to continually remesh after some iterations in order to produce a non-spurious mesh for the continuation of the calculation. Evidently, this re-meshing process can be heavily time consuming and complex thus the advantage of using an ALE mesh with respect to a Lagrangian one can be partially lost. where Finite Element Methods are used to discretize both the fluid and solid sub-domains and the aim is to simulate explosive charges contained within reservoirs with deformable walls. Within this study, the structure exhibits large displacements as well as material non-linearity. In order for the fluid ALE mesh to cope with the large displacements of the solid interface remeshing algorithms are employed to increase the accuracy of the coupling method.

[ VAN 07] carried out a comparative study of the ALE method applied to Finite Element discretizations with respect to the Fictive Domain method presented earlier. The results of this comparative study point to much more accurate simulations when the ALE method was used (in fact the ALE solution is considered the reference solution to which the Fictive Domain method solutions are compared). However, if active re-meshing is 1. Domain decomposition methods for fluid-structure interaction problems required, the ALE method loses some of its interest as the computational time is much longer than for the Fictive Domain method.

In [SOU 00], Souli et al. demonstrated the interest of using an ALE approach to fluid-structure interaction problems by showing academic and industrial examples where a Lagrangian approaches failed to converge due to excessive mesh distortions.

Lagrangian description for the fluid sub-domain

As already stated for the ALE mesh description, the use of a Lagrangian mesh-based method (FDM, FVM, FEM) would require complex numerical remeshing and meshrefinement techniques in order provide converged results. In the case of free-surface flows, the complexity of accurately tracking all the interfaces in question along with their possible fragmentations and reconnections would only complicate things further. For these reasons, a mesh-less based method for the spatial discretization of the fluid seems as the only plausible solution to overcome this challenge. By getting rid of the numerical mesh, the interfaces are tracked in a natural way as the material points or particles involved in the calculation spontaneously create the numerical support for spatial integration as the simulation goes along. Hence, no a-priori connectivity between particles is required and, in turn, the connectivity between the material points is determined at each new time-step.

Several types of mesh-free based methods exist [SUL 94, IDE 04, LUC 77, GIN 77], however this work is focused around a coupling done using specifically the Smoothed Particle Hydrodynamics method and the Finite-Element method. Hence, the SPH method will be the only meshless method presented within the scope of this work. The SPH method can also be used in and Eulerian and in an ALE framework as explained in [VIL 99] and as it will be discussed in Section 2.1. 3. However, considering the types of fluid flow and the dynamics of the structure in the simulations that will be featured in this work, only the expression of the SPH method in a Lagrangian framework will be of interest.

Concluding remarks on the different coupling approaches of semi-discrete numerical methods

During the present subsection, fluid-structure interaction coupling techniques have been presented according to different mesh descriptions for the fluid sub-domain (considering a Lagrangian mesh description for the solid sub-domain). It has been observed that when using an Eulerian description for the fluid, the tracking of boundaries can be quite complex if not impossible. Surface tracking can be made more approachable if an ALE mesh description is adopted, however if mesh distortions are too high, its efficiency can be hampered. This last statement is also applicable for mesh-based method expressed in a Lagrangian framework. Hence, it was put forward that only a mesh-less method can give the required versatility for carrying out the coupled simulations that will be carried out within the scope of the current work.

Domain decomposition methods -application to fluid-structure interaction problems

Techniques for breaking up large problems into collections of smaller and simpler to solve problems have been in use since more than a century [SCH 69]. The basic domain decomposition tools have been used in the past to couple problems which generally happen over the same physical continua and the decomposition helps apply different characteristics to different regions in the domain or to accelerate the solution process by superposing a collection of simpler broken down problems. Regarding fluid-structure interaction problems, these techniques can also be applied if the appropriate kinematic and dynamic conditions are applied at the interface as initiated by [BEL 80, COM 80]. The purpose of the current subsection is to present the basics behind the domain decomposition that inspired the fluid-structure interaction technique that is used within the scope of this work. Only the generalities of the domain decomposition methods will be given in the present section.

The specificities of how this technique is applied to couple a fluid and a solid domain, each being discretized with different techniques for space and time, will be explained in depth in the following chapter.

In general, when using domain decomposition techniques, the larger initial domain Ω is broken down into smaller sub-domains and Ω 2 separated by interface Γ 12 . On the right decomposition, the interface's dimension is the same as the dimension of Ω. On the decomposition in the left, the dimension of the interface is decreased by one [SCH 69] and are less commonly applied when solving coupling problems with semidiscritized methods since their implementation often calls for total control in the problem's matrix numbering in order to isolate the degrees of freedom (DOF) on the interface calculation points as a way to simplify the solution procedure of the resulting algebraic system [HER 02].

Ω i such that Ω = i Ω i .
Alternatively the Arlequin or Handshake method proposed by [DHI 98] can be used to couple problems in which different mechanical models describing different zones of the same domain are superimposed. Using this method, a duplication of mechanical states takes place along the overlapping area. Hence, in order to avoid duplicating the mechanical energy of the global system, the energy (virtual work associated to different forces) in the overlapped area is distributed using weight functions. Finally, in order to assure continuity of the mechanical state in the overlapping zone, bonding techniques such as the penalty method or the Lagrange multiplier method can be used.

Considering solid mechanics static problem split into two sub-domains each with an isotropic elastic material law, the Arlequin formulation allows to determine the displacement fields u 1 and u 2 that are kinematically admissible and the Lagrange multiplier λ allowing to minimise the functional in equation (1.4).

             W arl = W 1 (u 2 ) + W 2 (u 2 ) + C (λ, u 1 -u 2 ) W i (u i ) = 1 2 Ω i α i σ(u i ) : ε(u i ) dΩ - Ω i β i f i • u i dΩ ∀i ∈ (1, 2) C (λ, u) = 1 2 S λ • u + l 2 ε(λ) : ε(u) dΩ (1.4)
In equation 1.4,f represents the external forces, ε(u) is the strain tensor and σ(u) is the stress tensor associated to the displacement field u, α and β are the two weight functions (comprised between 0 and 1) associated to the internal and the external forces and l is a strictly positive parameter homogeneous to a length.

As stated before, this method becomes interesting when dealing with overlapping domains with different physical natures, such as the example put forward in [DHI 98] where a cantilever beam is split into a global domain, where typical beam theory can be used, and a smaller 2-D sub-domain for the region near the blocked boundary condition, where more localized mechanical effects can be studied.

From a fluid-structure interaction point of view, the Arlequin method seems difficult to put into practice since there is no overlapping zone for the fluid and solid sub-domains to share during a fluid-structure interaction problem. In the following sub-section, domain decomposition methods where no overlapping of sub-domains takes place will be presented and it will be seen how this type of coupling methods can be best suited to put into use when dealing with fluid structure interaction problems. Domain decomposition techniques in which no overlapping between sub-domains occurs can be simpler to implement and hence much more utilized in today's numerical simulations. Furthermore, this type of decomposition allows, in some cases, for problems in each sub-domain to be solved for independently, thus providing an ideal choice for parallel implementations. In the current section two types of decomposition techniques with no overlapping will be briefly presented : Schur's primal method and Schur's dual method. For a full run-down on other types of domain decomposition methods with no overlapping, the reader can consult: [DUR 97, CHA 14, FAR 91].

In order to facilitate the description of these two domain decomposition methods, the reference problem represented in Figure 1.11 will be considered.

From Figure 1.11, the domain Ω = Ω 1 ∪ Ω 2 ∪ Γ 12 . Each decomposed sub-domain has its own set of conditions acting on their respective boundaries: for i ∈ [1, 2], ∂ u Ω i for kinematic boundary conditions (Dirichlet-type) and ∂ F Ω i for dynamic boundary conditions (Neumann-type).

Domain decomposition methods for fluid-structure interaction problems

Ω 1 Ω 2 Γ 12 Ω 1 ∂ u Ω 1 ∂ F Ω 2 ∂ u Ω 2 ∂ F Figure 1.11: Domain decomposition reference problem
According to [HER 02, DUR 97], for the whole system Ω in Figure 1.11, considering a material with perfectly elastic behaviour, the variational form on displacement can be expressed as:

      ∀δU ∈ U O δ 1 2 Ω Tr[Kε(U)ε(U)]dΩ - Ω f d • UdΩ - ∂ F Ω F d • UdS = 0 U | ∂ u Ω = U d (1.5)
Where U is the unknown displacement vector and U O is the kinematically admissible field for the solution vector. K is Hooke's constant, since we are considering a perfectly elastic solid material for the sake of simplicity. For a static problem, finding the minima of (1.5) leads to the algebraic system of the following sort:

KU = F (1.6)
with K being the stiffness tensor, U containing the nodal kinematic unknowns and F representing the nodal applied forces on the system.

Remark: For the sake of simplicity in presenting domain decomposition methods, only the static analysis of mechanical problems will be considered in the current section despite the fact that the transient analysis implementation of these methods is the backbone of the fluid-structure interaction method used in this work. Further details regarding the implementation of the transient domain-decomposition methods will be covered in the next section and following chapters, where the fluid-structure interaction method will be covered.

Schur's primal method

Schur's primal method is based on establishing kinematic continuity at the interface between sub-domains. Hence, it's referred to as primal since it only makes use of the kinematic or primal variables to couple the different sub-domains involved. Dynamic equilibrium is verified later on resulting as a consequence of kinematic continuity conditions. Using the domain decomposition presented in Figure 1.11, the kinematic vector U can be considered a superposition of U 1 , U 2 and U Γ and by renumbering the DOF of the initial problem in a way that each of the sub-domain's internal DOF's are listed first followed by the interface's DOF's, (1.6) can be expressed as:

   K 11 0 K 1Γ 0 K 22 K 2Γ K Γ1 K Γ2 K ΓΓ       U 1 U 2 U Γ    =    F 1 F 2 F Γ    (1.7)
Using Gauss elimination's method, System (1.7) can be solved by solving first for the kinematic unknowns along the interface U Γ . This solution can then be found for the other two components of the unknown kinematic vector by substituting the value found for U Γ into the initial algebraic system (1.7). The full detail regarding the solution procedure of such a system can be found in [ DUR 97,CHA 14] Since this domain decomposition technique calls for the renumbering of the DOF in order to carry out the solution process of the algebraic system of equations (1.6), the implementation of this technique might call for an intrusive handling of the solvers used by the different sub-domains. Additionally for fluid-structure interaction problems, coupling relationships generally involve the relationship between the pressure at the interface and a kinematic variable (the velocity, in general). Thus a domain decomposition method in which the dynamic variables at the interface are also taken into account can be preferable. These latter type of decomposition techniques will be presented in the following sub-section.

Domain decomposition methods for fluid-structure interaction problems

Schur's dual method

The dual method for Schur's algorithm is one of the most widely used domain decomposition methods in solid mechanics [GRA 01, FAR 91, PRA 04]. As its name implies, this method makes use of a kinematic variable and dynamic variable (the kinematic's dual variable) in order to simultaneously enforce kinematic and dynamic equilibrium at the interface between sub-domains.

To establish Schur's dual method, once again, the initial domain Ω can be split into two sub-domains Ω 1 and Ω 2 separated by interface Γ 12 (Fig. 1.11). By following [FAR 91], one can build a Lagrangian functional based on (1.5) for each subdomain, J(v 1 ) and J(v 2 ) respectively, and by the imposure of a kinematic continuity condition at the interface written in its dual formulation with the help of a Lagrangian multiplier µ:

(v 1 -v 2 , µ) = Γ 12 µ • (v 1 -v 2 )dS (1.8)
Hence giving :

J * (v 1 , v 2 , µ) = J(v 1 ) + J(v2) + (v 1 -v 2 , µ) (1.9)
The solution of the minimization problem for the functional in (1.9) implies finding saddle point for J * (v 1 , v 2 , µ). It can be shown that the solution to this problem gives the two kinematic variables (u 1 and u 2 ) and the Lagrange multiplier (λ) guaranteeing kinematic and dynamic equilibrium at the interface.

As in the previous section, a static solid mechanics problem with a perfectly elastic constitutive law will be considered for the sake of simplicity. When using a semidiscretization technique, the solution to the minimization of the functional in (1.9) can be translated as the solution of the algebraic system (1.10). When this method is discritized using the Finite-Element Method, it is known as the Finite-Element Tearing and Interconnection method (FETI method).

       K 1 U 1 + L T 1 λ = F 1 K 2 U 2 + L T 2 λ = F 2 L 1 U 1 + L 2 U 2 = 0 (1.10)
If one compares this system with the problem stated in the previous section, the differences include the addition of the Lagrange multiplier vector (λ) as well as the constraint matrices (L i ) which can be considered as geometrical operators allow to project the DOF of each sub-domain unto the more restricted sub-space where the interface problem is solved.

In matrix terms, the system presented in (1.10) can be expressed as:

   K 1 0 L T 1 0 K 2 L T 2 L 1 L 2 0       U 1 U 2 λ    =    F 1 F 2 0    (1.11)
Following a bordered approach, the system in (1.11) can be solved by establishing the solution as the superposition to a constrained (a problem where the only applied forces on the domain come from the dynamic equilibrium at the interface) and unconstrained problem (a problem where only the external boundary and body forces on each separate sub-domain are considered). The details regarding the way system for this type of problem (1.11) is solved can be found in [FAR 91, GRA 01] with an extension to mismatching interface DOF's given by [HER 02]. The extension of this method to solve fluid-structure interaction method was developed in [ ZHE 13] and was extended to be used with different time scales in the current thesis.

In contrast to Schur's primal decomposition method, the method exposed in this subsection does not require an intrusive manipulation of the degrees of freedom on each sub-domain in order to provide an ideal way to find a solution for the algebraic system. Thanks to the constraint matrices (L i ), this method allows each sub-domain to be handled independently and only the continuity condition at the interface and the projection of the effect of the Lagrangian multipliers into each domain is made possible through the use of these geometrical operators. Furthermore this method ensures dynamic and kinematic continuity along the interface in a synchronous manner.

Multi-scale methods in time

The FETI methods, presented in the previous section, were introduced as powerful iterative solvers and remain one of the most used domain decomposition methods for static solid mechanics problems. These methods were later extended to dynamic problems in [FAR 95]. In the FETI method, a given spatial domain is divided into non-overlapping zones where an incomplete solution of the primal field (kinematic quantities) can be evaluated using a direct or iterative solver. Subsequently, field continuity is enforced via Lagrange multipliers applied at the interfaces of the concerned domains. This domain linking phase generates a symmetric dual problem where the unknowns are the Lagrange multi-

Regarding a widely used micro-time scale coupling technique, Gravouil and

Combescure proved in [GRA 01, COM 02] that imposing velocity continuity at the interface leads to a stable algorithm. In addition, they proposed a FETI-based multi-time-step coupling technique, known the GC method, which allows to couple arbitrary Newmark schemes with different time-steps in different subdomains. In this context, the GC method is shown to be unconditionally stable as long as each of the involved sub-domains satisfy their own stability requirements. However, for multi-time-step cases the GC method can induce energy dissipation at the interface, in contrast to the case where a unique time-step is used in all sub-domains, where the GC method is energy preserving. The details behind the GC method will not be developed in the current section since they will be presented in Chapter 3 of the current thesis.

Regarding a macro time-step based coupling algorithm, an interesting method has been recently proposed by Mahjoubi, Gravouil and Combescure called the MGC method [MAH 11]. The MGC method allows to couple sub-domain solvers with any time integrator even, for incompatible time-steps. This method is able to couple Simo, Krenk, Verlet, HHT and Newmark time schemes, thus, one of the great advantages of the MGC method is to be a general method able to couple a variety of time schemes. The MGC method holds a great advantage over the GC method which is that there is no energy dissipation at the interface and that it only requires one resolution of the condensed problem in a macro time step in order to obtain the Lagrange multipliers. The MGC method is a general but more complex method. It is solved not in velocity at a fixed time-step but in velocity increments and, particularly, the calculations of condensation operator H and of vector the right hand side operator b, deduced from the unconstrained quantities (not involving the influence of the sub-domain interface interactions) of both sub-domains are more difficult to obtain than for the GC technique. Furthermore, its implementation generally demands a more intrusive approach within the sub-domain specific solvers with respect to what can be achieved with the GC method. For this reason the MGC method was initially tested for implementation in the current work, but despite some early results, it was finally abandoned in favour of the GC technique whose micro time-scale coupling nature allows to better capture the non-linear effects of a rapidly changing fluid-structure interface.

One can add that the aforementioned transient dynamics coupling methods are endowed with advantages and drawbacks as, for example, the macro-based methods turned out to be more suited for parallel computing. However, since the interface constraint condition is only verified at the macro time-scale, an accurate tracking of the interface motion becomes more challenging. The micro-based multi-time-step coupling techniques handle more easily strong geometrical and material non-linearities at the fine time scale. However, some small amount of numerical dissipation can occur at the interface, due to the required interpolations of quantities at the large time scale. The present work targets the simulation of highly non-linear phenomena, both from the geometrical and the material perspectives, occurring in fluid-structure interaction problems. The GC method, due to its micro time-scale coupling characteristics, has thus been chosen in the current work and the details of its implementation will be presented and developed in Chapter 3 of the current thesis. The reader can consult [GRA 01, FAR 91, PRA 04, MAH 11] to obtain more information regarding these methods outside a fluid-structure interaction framework.

Conclusion on the domain decomposition methods and its applicability to fluid-structure interaction problems

The domain decomposition methods presented in the current subsection were originally intended for their application in solid mechanics problems [DUR 97, FAR 91]. As previously discussed, these decomposition techniques allow to simplify the solution process of a statics problem in which different constitutive models or boundary conditions can be used in each sub-domain making part of the original problem. In transient solid dynamics, this type of decomposition techniques can be shown to be applicable [GRA 01, PRA 04] and their interest can be magnified if the different sub-domains that make-up the original domain require different time-scales in order to capture different time depending characteristics that each sub-domain may present.

Domain decomposition methods for fluid-structure interaction problems

Hence, the extension of the domain decomposition techniques onto transient fluidstructure interaction problems can be envisaged if the whole fluid-solid system is considered as an initial domain that is decomposed into a fluid and solid sub-domains. The dual-Schur method offers greater flexibility within the large family of domain decomposition techniques, as there is no need to interfere with the way in which the DOF of each sub-domain are organized since especially built constraint matrices allow to extract the fluid-solid boundary DOF from each sub-domain in order to solve the continuity equations at the interface and to project back the solution onto each sub-domain. Furthermore, the dual-Schur method, in contrast to the primal method, gives an insight into the forces that operate at the interface linking both domains as well as other involved kinematic quantities. Since in many cases boundary conditions at the fluid-solid interface for transient fluid mechanics problems involve a relationship between the pressure and the velocity at the interface, the fact that both dynamic and kinematic variables are present in the dual-Schur method provides an entry point the coupling that can be developed for both sub-domains [ ZHE 13]. As already exposed in the current subsection, all the details regarding the implementation of the dual-Schur domain decomposition method will be presented during the following sections of the current work. An extension on applying the method to problems featuring different time-scales in the solid and the fluid with the presence of strong non-linearities will also be given.

Solid contact-impact modelling

Numerical modelling of impacts between solids is quite a vast research subject and several techniques exist which attempt to recreate as accurately and inexpensively as possible the complexity of this recurrent physical phenomenon [HAL 85, BEL 91]. For the current work, an appropriate way of simulating solid impacts is necessary in order to be able to predict the extent of the damage that can be caused on a functioning hydroelectric turbine blade by an unexpected stone found within the jet stream. Besides being accurate and reliable, this method must be made compatible with the fluid-structure interaction simulation strategy that is used to link the fluid and solid sub-domains. The numerical methods for simulating impacts using the finite-element method will be handled by a specialized solver and the basic characteristics behind the algorithm taking care of this complex task will be briefly presented in the current subsection. The implementation algorithm making the impact modelling compatible with the sub-domain coupling technique will be given Solid contact-impact modelling in the following sections.

Impenetrability conditions

For the sake of simplicity two solid bodies, Ω 1 and Ω 2 can be considered which can be discterized using a Lagrangian mesh (mesh point motion is governed by the continuum's governing equations). In order to guarantee that there will be no interpenetration of the two bodies, it must be ensured that neither body can occupy the same space at the same time.

Ω 1 ∩ Ω 2 = / 0 (1.12)
Although it is impossible for Ω 1 and Ω 2 to occupy the same physical position at the same time, the two domains can be in contact over a surface Γ c as Figure 1.12 shows. Over surface Γ c , a mathematical condition to ensure impenetrability between the nodes and elements of each domain must be established. In general, a condition involving the normal velocities to the interface as well as a dynamic condition involving the normal traction to the surface are considered. Conditions on the local tangential components of the velocities and tractions can also be established in the case in which friction is considered for the contact. However, for the sake of simplicity, we will only consider the local normal components for this conditions which read : 1. Domain decomposition methods for fluid-structure interaction problems

γ n = (v 1 -v 2 ) • n 1 ≤ 0 t n = t 1 N (x,t) = -t 2 N (x,t) ≤ 0 (1.13)
The first condition in equation (1.13) guarantees velocity equilibrium along the normal to body Ω 1 on contact domain Γ c . The fact that this equation is expressed as an inequality, emphasizes that when both bodies are in contact, they can either remain in contact (γ n = 0) or separate (γ n < 0). For the second condition, t i N (x,t), i ∈ [1, 2] are the normal tractions expressed on time t and on a point in space x lying on Γ c . The normal tractions are obtained through projections of the Cauchy stress tensor on the normals to the interface. This latter condition is known as the traction equilibrium condition and it only takes into account compressive tractions (adhesion between surfaces is not considered for simplicity). Both of these conditions can be combined in order to express a unitary contact condition:

γ n • t n = 0 (1.14)
Condition (1.14) expresses the fact that when the two bodies are in contact, they must either remain in contact (γ n = 0, t n < 0) or separate (γ n < 0, t n = 0).

Contact detection and contact force determination using the FEM

When using a Finite-Element method to simulate impacts between solids, several exist techniques which are capable of detecting zones of contacts and subsequently providing the contact forces acting on the solids and preventing interpenetration. Two of the most widely employed algorithms used today are the sliding lines/surfaces algorithm and the pinball algorithm [HAL 85, BEL 91].

The sliding lines/surfaces algorithm (Figure 1. [START_REF] Zhe Zhe L | Développement d'une méthode de simulation de couplage fluide-structure à l'aide de la méthode SPH[END_REF]) is one of the first ones to be used for simulating impacts between solids using the FEM method among others. It directly uses the mesh elements and nodes in order to carry out geometrical operations allowing to determine the existence of contact between solid bodies discretized in space using a Lagrangian mesh. The method is based on a "master" and "slave" pre-determination of bodies under potential contact in which the nodes of the master surface are tested for proximity with nodes of the "slave" surface at each time iteration. Interpenetration is detected by checking each "slave" node's position against the surface or edge around each "master" node. If interpenetration is detected, resulting forces can be calculated using several methods such as the penalty method, or using Lagrangian multipliers. The full details regarding the sliding lines algorithm can be found in [HAL 85]. The pinball algorithm (Figure 1.14) was developed in order to provide a solution to some of the shortcomings of the sliding lines/surfaces algorithm such as the complex geometrical calculations it requires in order to detect contacts between domains or its inability to detect certain types of contacts between edges and surfaces of the involved meshes. This more recently developed algorithm relies on a pseudo re-meshing of the surface elements in the contacting bodies with circles/spheres whose sole purpose is to detect a contact through simple mathematical condition regarding their radii. The possibility to refine the density of the pinballs without needing to change the original mesh is one of its main advantages when compared to the sliding lines/surfaces algorithm. The calculation of efforts at the interface following a contact between bodies can also be done through the penalty method or using Lagrange multipliers. The full details regarding the formulation and implementation of the pinball method can be found in [BEL 13, CAS 02].

Following contact detection using the techniques described in the previous section, the contact force impeding solid interpenetration needs to be obtained in order to obtain the external efforts acting on each solid and thus determine their new configuration. As put forward in [HAL 85, BEL 91, CAS 02], several methods to determine contact forces exist such as the penalty method and the Lagrange multiplier method. For the current work, we will only focus on giving the details on the effort determination through the use of Lagrange multipliers as it is the same method that is used for determining the fluid-induced efforts on the structure.

For the current work, the use of precise algorithms for simulating contacts between solid bodies is necessary. The solid solver, Europlexus, that will be used for the simulations in this work, which is based on the Finite-Element method provides two algorithms capable of calculating the contact forces between solids. These algorithms are known as the sliding lines algorithm and the pinball algorithm. The calculations of solid contacts will be considered as part of a single solid calculation even in the presence of more than one solid domain. Hence the external forces produced through contact will be part of an monolithic calculation within the solid solver and the rest of the external efforts on the solid sub-domains will come from the fluid-structure interaction through a partitioned interaction with an external SPH-based fluid solver. Such an approach involving the decoupling of solid-solid calculations and the use of domain decomposition methods can be seen in the works of [BOU 06, GRA 01] and provides an elegant solution to coupling external solvers while taking into consideration impacts between solid sub-domains.

Conclusion

Over the course of the current chapter, several techniques to solve fluid-structure interaction problems have been presented and their relevance with regards to the objectives of the current work has been evaluated. Domain decomposition methods allowing to split a mechanical problem into separate and simpler to solve entities have also been presented.

As it has been explained, these decomposition techniques, although initially conceived to be used in solid mechanics statics problems, can be extended to be used in fluid-structure interaction problems if considering the fluid and solid as separate sub-domains of a whole initial problem. Finally, numerical methods to simulate contacts between solids have also been introduced as their understanding is important in order to carry out the simulations within the scope of this work. From all the information presented, the most convenient way of carrying out the simulations has been argued to be a coupling between the Finite Element Method for the simulation of the solid and the SPH method for the fluid sub-domain using a domain decomposition method in which both sub-domains are linked using a Dual-Schur-inspired (FETI) method. Solid contact simulations will be handled as internal calculations within the solid solver.

The following chapters in this text will focus on precisely describing the governing equations for each sub-domains along with the basic concepts around the semidiscretization techniques used to solve them. The coupling technique will be described in detail and validation examples will be given. Subsequently, using the numerical coupling method described, calculations of the damage induced by a rock impact on a hydraulic turbine blade will be carried out in order to show the robustness of the proposed method.

Chapter 2

Fluid-solid coupling strategy

In this chapter, the fluid and solid sub-domain governing equations will be presented along with their respective semi-discretized forms providing the numerical framework for the numerical simulations done within the scope of this work. Subsequently the domain decomposition technique allowing to couple both discritized sub-domains will be presented In the current section, the equations governing the fluid sub-domain are presented first along with the numerical mesh-less technique used to discretize the physical domain. Special attention should be paid to the boundary condition treatment, which allows to set up a system of coupled equations whose unknowns are the pressures and velocities at the interface with the solid sub-domain.

Governing equations for the fluid

For the fluid sub-domain, water is the only liquid considered. Hence an equation model capable of describing the behaviour of a Newtonian, weakly compressible fluid such as the Navier-Stokes model [TEM 84] is to be applied. Furthermore, only isothermal flows at high Reynolds numbers will be considered (inertial factors are thus considered preponderant over viscous factors), giving way for the energy conservation equations as well as all viscosity-related terms to be discarded from the Navier-Stokes model. In this manner, the Euler model of the conservation equations is to be considered (2.1) for the description of the behaviour of fluid flows in the current study.

         ∂Φ ∂t + ∇ • F c = Q Φ = ρ f ρ f v f , F c = ρ f v f ρ f v f ⊗ v f + p f I , Q = 0 ρ f g (2.1)
In equation (2.1), the Euler model of the conservation equations are expressed using the vector Φ which is known as the vector of conservative variables ρ f , the fluid density, and ρ f v f , the momentum , which is the product of the fluid's density and the fluid's velocity v f . F c is called the convective flux tensor and it can be expressed in terms of Φ if considering a weakly compressible state equation allowing to express the fluid pressure p f (I being the identity matrix) in terms of the fluid density ρ f . Finally, the vector Q contains the source terms of the conservation equations, which in the current case, only consider the gravitational induced body forces through the gravitational field from which the gravity vector g is obtained.

The differential form of the Euler equations for fluid flows is valid for functions that are at least C 1 since the divergence of the convective flux tensor is needed to solve the system of equations. For cases in which the conservative variables are piecewise differentiable, it 2. Fluid-solid coupling strategy is interesting to express system in the integral form over an infinitesimal fluid volume Ω f and using the Reynolds Transport Theorem as presented in (2.2)

∂ ∂t Ω f Φ dΩ f + Ω f ∇ • F c dΩ f = Ω f Q dΩ f (2.2) Expression (2.
2) is explicitly expressed in integral form as :

         ∂ ∂t Ω f ρ f + ∇ • ρ f v f dΩ f = 0 ∂ ∂t Ω f ρ f v f dΩ f + Ω f ∇ • ρ f v f ⊗ v f dΩ f = Ω f p f IdΩ f (2.3)
The unknowns in (2.3) are ρ f , p f and v f . The third equation needed to close the system comes from the weakly compressible equation of state giving the relationship between the fluid pressure p f and the fluid density ρ f :

p f = ρ 0 f (c 0 f ) 2 γ f ρ f ρ 0 f γ f -1 (2.4) Equation (2.4
) is known as the Tait equation of state [MAC 66], where γ f ≈ 7 and the 0-superscripted terms are employed to denote the reference values of the pressure p 0 f the density ρ 0 f and the sound speed c 0 f .

The equations in (2.2) are expressed in an Eulerian framework since the reference volume Ω f is fixed in space and the fluid material points move through this volume serving as support for the evaluation of the integrals. In order to consider the concept of the motion of the reference volume a new parameter v 0 should be introduced allowing to obtain the material time derivative of the points within the moving reference volume. In such a case [HIR 74, NEU 14], the conservation equations can be expressed as :

∂Φ ∂t + ∇ • (v 0 ⊗ Φ) dΦ dt +∇ • (F c -v 0 ⊗ Φ) F(Φ,v 0 ) = Q (2.5)
(2.5) can be expressed in integral form over the Lagrangian reference volume, Ω L f , moving with velocity v 0 as:

r ij
d dt Ω L f Φ dΩ L f + Ω L f ∇ • F(Φ, v 0 ) dΩ L f = Ω L f Q dΩ L f (2.6)
In expressions (2.5) and (2.6), the velocity of the reference volume v 0 can take different values in order to find the expression of the governing equations under different configurations. If v 0 is set equal to zero, the governing equations are expressed in an Eulerian framework. If v 0 is set to be equal to fluid flow velocity, a Lagrangian configuration of the governing equations is obtained. Finally, a third configuration in which, the velocity of the reference volume v 0 is set arbitrarily in order to suit some constraint inherent to the calculation under consideration, allows to express the conservation equations in an arbitrary Lagrangian-Eulerian (ALE) framework.

Discretization using the SPH method

The equations presented in (2.3) will be solved using a mesh-less technique known as Smoothed Particle Hydrodynamics. This technique is based on the use of a compactly supported Gaussian-like function W (r, h), called the kernel function, used to weigh interactions between the unconnected numerical particles discretizing a physical domain.

Figure 2.1 shows the shape of a typical Gaussian-like kernel function used in SPH. The kernel function W (r, h) is centered on a numerical particle i, where r i j represents the distance separating it from another arbitrarily positioned particle j and h is a parameter used to find the radius of the function's support domain D i bounded by ∂D i . Outside the support domain of the function, its value is nil. For the current study, we will make use of the Wedland C4 kernel function [SCH 06]. There are several families of SPH-based techniques that can be employed to simulate different types of fluid flows.

The traditional formulation of the Smoothed Particles Hydrodynamics method as proposed by [LUC 77, GIN 77] is written in a Lagrangian framework in a way guaranteeing at least zero order consistency for the spatial approximation of the physical domain under consideration and using an artificial viscosity term Π i j to stabilize the numerical scheme.

Using the original SPH formulation, the semi-discrete fluid conservation equations give :

                   dx i dt = v i dρ i dt = ρ i ∑ j∈D i (v j -v i ) • ∇ j W i j ω j dv i dt = g i -∑ j∈D i m j p j ρ 2 j + p i ρ 2 i + Π i j ∇ j W i j
(2.7)

In system (2.7), an equation has been added to the conservation equations in order to keep track of the position x i of a numerical particle i. The discretized conservation equations, consider the interaction between a particle i and its neighbour j under the support domain of D i . For simplicity, equations in system (2.7) only take into account particle interactions far away from the boundaries of the discretized domain and no truncation of the support of the kernel function ∂D i is considered. Notice that using the SPH method, the calculation of spatial gradients is done on the weight functions and not directly on the conservative variables mimicking the purpose of the shape functions used in the finite-element method.

SPH-ALE method

Despite its ease of implementation and proven accuracy for certain cases [GIN 77, DEL 11], the traditional Lagrangian SPH method will not be the one used to discretize the fluid sub-domain in the current work. Instead, a derived ALE formulation proposed by [VIL 99] will be used since it provides major flexibility as well as a numerical framework for elegantly imposing boundary conditions which make the SPH-ALE method the proper choice for dealing with complex solid geometries [MAR 07].

Using the notations in (2.5) and the divergence theorem, the conservation equations discretized using the SPH-ALE method are expressed in semi-discrete form for any particle i within the computational domain as :

d(ω i Φ i ) dt + ω i ∑ j∈D i ω j (F i + F j ) • ∇ i W i j + ω i ∑ k∈∂D i s k (F i + F k ) • n k W ik = ω i Q i (2.8)
where F i and F j represent the convective fluxes evaluated at the positions of particles i, x i , and j, x j , respectively, as well as the convective flux tensor of boundary element k of surface s k and at position x k .

The sum of convective fluxes on different particles in equation (2.8) guarantees conservation of momentum between particles since the influence of particle i on particle j can be shown to be the opposite of the influence of particle j on particle i [DEL 11]. Nevertheless the formulation in (2.8) doesn't ensure mass conservation between particles at a local scale even if mass is conserved at the global scale. Furthermore, as seen in [VIL 99] the semi-discretization of the conservation equations in (2.8) is unconditionally unstable. An artificial viscosity as shown in the original Lagrangian formulation of the SPH method could be a possibility for providing stability to the numerical method. Nevertheless, another approach was presented by Vila in [VIL 99], in which the summation of fluxes between particles, F i + F j , is replaced by intermediate state between them, 2G E,i j = 2G E (Φ i , Φ j ), which is calculated as the solution to a Riemann problem between the states of particles i and j. Using this new formulation, equation (2.8) can be expressed as :

d(ω i Φ i ) dt + ω i ∑ j∈D i ω j (2G E,i j ) • ∇ i W i j + = ω i Q i (2.9)
where:

G E,i j = ρ E i j (v E,i j -v 0,i j ) ρ E i j v E,i j ⊗ (v E,i j -v 0,i j ) + p E,i j I (2.10)
and where the wall terms k ∈ ∂D i from (2.8) are excluded for the sake of simplicity.

In equation, (2.10), the top line is the numerical mass flux expressed from the solution of the Riemann problem for the mass conservation equation. The bottom line is the numerical flux expression to the solution of the Riemann problem in the case of the momentum conservation equation. The former and latter expressions are in turn composed of 2. Fluid-solid coupling strategy the solution to the Riemann problem to the density ρ E i j , the velocity v E,i j and the pressure p E i j whose solution process will be explained in the following lines. One should notice the arbitrary velocity v 0,i j = v 0 (x i j ,t) with x i j = (x i + x j )/2 which takes into account the movement of the numerical particles independently from the flow velocity of the fluid continuum (ALE-nature of the formulation).

In order to obtain the solution to the intermediate state between two particles i and j, a local axis system is considered for a 2D problem in which the 1st direction is tangential to the vector pointing from i towards j and the 2nd direction is orthogonal creating a 2D direct local orthonormal system (O, x ′ , y ′ ) as portrayed in Figure 2.2. By establishing the new coordinate system, the solution to the conservation equations will only be sought on the direction parallel to x ′ at half the distance between the x i and x j by solving a Riemann problem [TOR 13] at the numerical discontinuity between the two disjoint numerical entities. The contribution along the tangential direction y ′ will be obtained when solving the Riemann problem between particle i and other surrounding particles different than j. Thus, the conservation equations are projected on a 1D vector subspace simplifying the solution process.

Following [MAR 07, ZHE 13, TOR 13], after some mathematical simplifications, the algebraic system that must be solved for in order to find the value of the primitive quantities (ρ, v x ′ , v y ′ ) when the problem is projected on x ′ is :

             ∂W ∂t + A ∂W ∂x ′ = Q W =    ρ v x ′ v y ′    , A =    v x ′ ρ 0 c 2 /ρ v x ′ 0 0 0 v x ′    , Q =    0 g x ′ g y ′    (2.11)
In equation (2.11), c represents the speed of sound in the fluid, v x ′ and v y ′ represent the components of the velocity vector v in directions x ′ and y ′ , respectively. Similarly, for the gravity vector g, g x ′ and g y ′ represent the components in their respective directions.

The solution for the non-linear hyperbolic PDE system in (2.11) is obtained by solving a Riemann problem whose solution is a function of the known states of the primitive variables in particle i and particle j. The pressure p variable, is obtained through the state equation (2.4). To set up the Riemann problem between particle i and its neighbour j, the characteristic waves associated to the eigenvalues of system (2.11) are needed.

These are associated to the eigenvalues of matrix A, hence, the roots of the characteristic polynomial det(A -λI) = 0 are retrieved:

       λ 1 = v x ′ -c λ 2 = v x ′ λ 3 = v x ′ + c
(2.12)

The eigenvalues in (2.12) represent the speed of propagation of information given and initial value in time and space. As presented in Figure 2.3, one can plot the eigenvalue in a x ′ -t plane in order to see the different domains on which the sought solutions to the Riemann problem may exist.

= ρ * R = ρ * and v * L,x ′ = v * R,x ′ = v * x ′ ).
The remaining two waves associated to eigenvalues λ 1 and λ 3 can either be a shock-wave or a rarefaction wave. The values of the conservation variables are known for W L and W R since they represent the states of particles i and j or a linearization of these states. The value of the conservative variables needs to be obtained for the * zone in order to determine the nature of the waves associated to the remaining eigenvalues. An algebraic relation for the Euler equations, matrix A contains non-constant values and hence the hyperbolic system of PDE in (2.11) is non-linear. In order to simplify the solution process, the matrix A in system (2.11) can be linearized as follows:

A =    v x ′ ρ 0 c 2 /ρ v x ′ 0 0 0 v x ′    , where              ρ = 1 2 (ρ L + ρ R ) v x ′ = 1 2 (v L,x ′ + v R,x ′ ) c = 1 2 (c L + c R ) (2.13)
Using the approximation in (2. [START_REF] Zhe Zhe L | Développement d'une méthode de simulation de couplage fluide-structure à l'aide de la méthode SPH[END_REF]), solving the approximate Riemann problem provides the following relations between the known states of the conservative variables and the states of these variables in the * region:

       ρ * = 1 2 (ρ L + ρ R ) + ρ 2c (v L,x ′ -v R,x ′ ) v * x ′ = 1 2 (v L,x ′ + v R,x ′ ) + c 2ρ (ρ L -ρ R ) (2.14)
Furthermore, since the tangential velocity component (v y ′ ) does not change in through the 1-D shockwave or a rarefaction wave associated to the eigenvalues λ 1 and λ 3 , the values are conserved between the initial domains and the * domain, hence :

v * L,y ′ = v L,y ′ v * R,y ′ = v R,y ′ (2.15)
The approximation allowing to find the expressions of the unknown variables in the star-region is called the Primitive Variable Riemann Solver (PVRS) and for the applications in the current work it has been shown to provide sufficient accuracy as well as a considerable time gain since there is no need to solve a non-linear problem between each pair of particles within the considered fluid domain [MAR 07, NEU 14].

After obtaining the values over the whole x ′ -t domain, the nature of the waves associated to eigenvalues λ 1 and λ 3 can be determined by comparing the density in the * region to the densities in the initial left and right states. Following this step, the final solution to the Riemann problem can be found by comparing the projected velocity of the reference volume v 0 along the unitary vector directed form particle i to particle j to the value of the eigenvalues inherent to the linearised Riemann problem between both particles. In the current work, the type of flows that will be considered will take into account low Mach numbers, hence the flow velocity will be considered much lower than the sound speed. For this reason, the solution to the Riemann problem will be considered to be found exclusively within the * region of the x ′t domain [ZHE 13, NEU 14, REN 15].

As such, the solution to the Riemann problem can take the following form :

       ρ E,i j = ρ * v E,i j = v * x ′ x ′ p E,i j = p * (2.16)
Furthermore, the initial states Φ L and Φ R can be calculated in different ways. The simplest and one of the earliest ways proposed to obtain the left and right states for the Riemann problem is to consider the values of the conservative variables over the whole domain up to the discontinuity x i j = (x i + x j )/2 as being piece-wise constant. Hence, the initial left and right values for the Riemann problem are automatically obtained and the solution of the problem is straightforward through (2.14) and (2.16). Despite the simplicity in using this approach proposed by Godunov in [GOD 59], the influence of the constant value field hypothesis can be rather strong and, as seen in [MAR 07, LED 10], can add considerable numerical viscosity to the numerical scheme, thus leading to spurious results. An alternative to overcome the drawbacks of the constant value field 2. Fluid-solid coupling strategy hypothesis, the Monotone Upstream-centered Scheme for Conservation Laws (MUSCL) technique proposed by [ VAN 79] suggests the reconstruction of the fields in the left and right initial states by using a linear interpolation :

Φ L = Φ i + α(Φ i , Φ j , ∇ i Φ)∇ i Φ • x i j Φ R = Φ j + α(Φ i , Φ j , ∇ j Φ)∇ j Φ • x i j (2.17)
where α(Φ i , Φ j , ∇ i Φ) represents a non-linear limiter function added to the linearised approximation in order to prevent spurious oscillations and ∇ i and ∇ j represent the gradients evaluated at x i and x j , respectively through the classical SPH approximations [GIN 77, VIL 99]. Taking into account the wall terms, the full set of governing equations semi-discretized through the SPH-ALE method becomes:

                                                     dx i dt = v 0,i dω i dt = ω i ∑ k∈∂D i W ik (v 0,k -v 0,i ) • n k s k + ω i ∑ i∈D i ∇ i W i j • (v 0, j -v 0,i ) ω j d(ρ i ω i ) dt = -ω i ∑ k∈∂D i W ik 2ρ E,ik (v E,ik -v 0,ik ) • n k s k -ω i ∑ j∈D i 2ρ E,i j (v E,i j -v 0,i j ) • ∇ i W i j ω j d(ρ i ω i v i ) dt = ρ i gω i -ω i ∑ k∈∂D i W ik 2[ρ E,il v E,ik ⊗ (v E,ik -v 0,ik ) + p E,ik I] • n k s k -ω i ∑ j∈D i 2[ρ E,i j v E,i j ⊗ (v E,i j -v 0,i j ) + p E,i j I] • ∇ i W i j ω j (2.18)
In (2.8), the arbitrary particle velocity vectors of particles i and its neighbour j, as well as wall element k are represented by v 0,i ,v 0, j and v 0,k , respectively.

Although the SPH-ALE technique allows to impose a numerical particle velocity that can be different from the fluid material velocity, we will use its Lagrangian version because we are interested in having an accurate tracking of a constantly evolving fluid-solid interface. Additionally, Vila's formulation of the SPH equations provides a necessary mathematical formalism that was subsequently exploited by Marongiu in order to find the value of the pressure at the solid wall by solving a linearised partial Riemann problem [MAR 07]. In order to obtain system (2.18) under a Lagrangian formulation, the arbitrary particle velocity vectors are set equal to the fluid velocity (v 0,i = v i and v 0, j = v j ). Since the velocity of the wall is always independent from the velocity of the fluid continuum, the velocity wall terms ins 2.18 can be expressed as v 0,k = v k . Furthermore, we will consider exclusively the case where the boundary of the fluid domain is an impermeable wall. We will consider v E,ik = v 0,ik , thus no mass transfer exists between the fluid sub-domain and the walls bounding it. Finally, we impose v i j = (v i + v j )/2, which implies that the interface between two particles where the solution to the Riemann problem will be sought, is imposed halfway between them. Let us observe that this point is of ALE nature because its velocity is not the same of that of the corresponding material point. Hence, the semi-discrete SPH-ALE system of equations expressed in a Lagrangian framework gives:

                                     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 i j (v j -v i )ω j d(ρ i ω i ) dt = -ω i ∑ j∈D i 2ρ E,i j (v E,i j -v i j ) • ∇ i W i j ω j d(ρ i ω i v i ) dt = ρ i g i ω i -ω i ∑ k∈∂D i W ik 2p E,ik • n k s k -ω i ∑ j∈D i 2[ρ E,i j v E,i j ⊗ (v E,i j -v i j )p E,i j I] • ∇ i W i j ω j (2.19)
Where x i is the position of a fluid particle i and ω i is its volume. The velocity vectors of particles i and its neighbour j, as well as wall element k are represented by v i ,v j and v k , respectively. Additionally, ρ i represents the density of particle i and p i represents its pressure. In this formulation, the boundary terms for a particle i, ∂D i , are taken into account as they will be subsequently used to properly treat boundary conditions.

In system (2.19), as it was previously stated, the terms ρ E,i j , v E,i j and p E,i j represent an intermediate status between two fluid particles that is obtained by solving a linearised Riemann problem at the interface between the particles. System (2.19) is more convenient for solving FSI problems since in a Lagrangian formulation, the interfaces -between the fluid and any other media -are tracked in a natural way.

Fluid-solid coupling strategy

Wall boundary conditions

The truncation of the support domain of the kernel function has always been a delicate matter when seeking to properly impose boundary conditions on a finite-domain SPH To achieve this, it was established that the term p E,ik in (2.19) represents the pressure solution of a partial Riemann problem that can be expressed between a fluid particle i and the boundary element k in its neighbourhood. When a solid boundary condition is imposed on the boundary element, the velocity solution v E,ik equals the physical boundary condition (usually the velocity of the solid) and p E,ik results from the combination of this imposed condition and the fluid state (pressure and velocity) of particle i. The method has proven to be efficient at imposing boundary conditions even when confronted with complicated solid geometries.

In order to obtain the expression of the value of the pressure on a solid wall element k, the boundary term of momentum equation must be considered (2.19-d). Let us call it b k :

b k = ∑ k∈∂D i ω i W ik 2p E,ik • n k s k (2.20)
where we consider, for a solid wall element k, all the fluid particles i whose support domain is truncated by k.

By writing that the opposite of the contribution to the change of momentum brought by the wall element k to the fluid (2.20) is equal to the hydrodynamic force exerted by fluid on element k we obtain:

p k = ∑ i∈∂D k ω i W ik 2p E,ik (2.21)
where the normal and surface terms n k s k have disappeared since they are present on both sides of the equation.

Using a linearised Riemann solver [TOR 13, DUB 01], we can find the solution for the pressure p E,ik by considering a fluid particle i and a solid element k. The quantity p E,ik is expressed as a function of the wall velocity v k :

p E,ik = p i -ρ i c i (v k -v i ) • n k (2.22)
Fluid boundary equations for fluid-structure interaction

Using the information given in the previous sub-section, the relation between the pressure at the wall and the velocity of the wall can be established. This expression is obtained by substitting (2.22) into expression (2.21) :

∑ i∈D k 2ω i ρ i c i W ik v f k + p k = ∑ i∈D k 2ω i (p i + ρ i c i v i • n k )W ik (2.23)
where v f k corresponds to the value of the solid element velocity in the outwardpointing normal direction, v f k = v k • n k . p k represents the value of the fluid pressure on the solid wall-element.

Finally, equation (2.23) can be written in matrix form as:

M f v f k + Λ = G f (2.24)
where M f denotes a square diagonal matrix of dimension (N k × N k ), N k being the total number of solid wall elements. For k ∈ (1, N k ), the coefficients along the diagonal are:

M f (k, k) = ∑ i∈D k 2ω i ρ i c i W ik (2.25)
the vectors v f k , Λ and G f are of dimension (N k × 1) and their coefficients for k ∈

(1, N k ) are given below:

         v f k (k) = v f k = v k • n k Λ(k) = p k G f (k) = ∑ i∈D k 2ω i (p i + ρ i c i v i • n k )W ik (2.26)
In order to integrate in time the fluid status, the 2nd-order accurate Runge-Kutta 2 midpoint scheme will be employed. As this integration scheme requires the determination of the fluid status at an intermediate stage (t n+1/2 ) to predict the corrected solution at the final stage (t n+1 ), the system (2.24) must be expressed at both time stages:

   M n+1/2 f v n+1/2 f k + Λ n+1/2 = G n+1/2 f M n+1 f v n+1 f k + Λ n+1 = G n+1 f (2.27)
Equation ( 2.27) provides the relationship between the solid boundary velocity and the pressure at the interface corresponding to the unknowns of the system for the fluid subdomain. The complete system of equations needed for the coupling between the fluid and solid sub-domains will be completed with the solid and interface equations that will be presented in the following sections.

Governing equations for the solid sub-domain

The reference problem from which the governing equations of the solid sub-domain can be derived is shown in Figure 2. 4. In the latter, a solid domain Ω whose frontier ∂Ω is subject to kinematic boundary conditions on ∂ u Ω and to dynamic boundary conditions on ∂ F Ω is depicted. Ω is subject to a body force induced by a gravitational force field.

Using the principle of virtual power [ZIE 77, BEL 13] and basic principles from continuum mechanics, one can establish the governing equations for the solid sub-domain depicted in Figure 2. 4. To simplify this derivation, small displacements are considered, meaning that the different solid configurations of the solid sub-domain are merged with the initial configuration. The unknowns of the problem are the displacement u(X ,t) and the stress tensor σ(X ,t) on every point of the domain, ∀X ∈ Ω at any t in time.

In order to solve the reference problem, one must find u(X ,t) and σ(X ,t) verifying in each point and at every instant the following:

                   ∀t, u(t) | ∂ u Ω = U d ∀X ∈ Ω, u(X ) | t=0 = U 0 , du dt (X ) | t=0 = V 0 ∀t, ∀ û ∈ Û, Ω ([ f d + ρ d 2 u dt 2 ] • û -Tr[σε( û)]) dΩ = ∂ F Ω F d • û dS ∀t, σ = Kε(u) (2.28)
where U d is the imposed kinematic boundary condition on ∂ u Ω, U 0 and V 0 are the initial displacement and velocity conditions respectively, û is the displacement test function that is kinematically admissible on Û, which is the subspace of test functions which impose a zero displacement on ∂ u Ω, ε(u) is the strain symmetrical tensor, σ is the stress symmetrical tensor, K is the stiffness matrix obtained for a linear elastic material in small displacements, ρ is the material's density, f d are the body forces the gravity field applies on the volume and F d are the dynamic conditions applied on ∂ F Ω.

Finite-Element method for the solid

In order to find a set of exact solutions for the conditions and equations presented in (2.28) it would be necessary to analytically solve the system in question. However, analytical solutions for hyperbolic problems such as the one presented previously rarely exist in practice. One can attempt to find an approximate solution of system (2.28) by discretizing the solid domain Ω into a domain of finite dimension Ω h where the solution u(X ,t) can be approximated by u h (X ,t). This is the cornerstone idea behind the Finite Element method which, as already announced several times, will be used for the simulation of the solid domain in the current work. Only a brief overview of the main ideas behind the Finite-Element method will be given in this work. For a full run-down of all the principles and hypotheses that are behind this method the reader can consult [ZIE 77, BEL 13].

The approximate solution u h (X ,t) can be expressed as the the sum over the n nodes used to build the finite elements discretizing Ω h of the product of a shape function φ and the discrete values of u at each of these nodes:

u h (X ,t) = n ∑ i=1 φ i (X )u i (t)
(2.29)

Fluid-solid coupling strategy

Using this discretization, the solution vector u h is a member of a vector subspace of finite dimension n () and thus an algebraic system can be constructed in order to attempt to solve system (2.28) in a approximative way. The system to solve for reads:

           ∀t, u h (t) | ∂ h u Ω = U h d u h | t=0 = U h 0 , du h dt | t=0 = V h 0 ∀t, ∀ ûh ∈ Û h , T ûh • (M s üh + F int -F ext ) = 0 (2.

30)

Where the h-superscripted variables and domains represent the discretized counterparts of the variables presented in (2.28) taking into consideration the restriction shown in (2.29). The third line of equation, presents the Mass matrix M s whose components on a i-th line and j-th column (for any given finite element e with volume V e ) are given by:

M i j = ∑ e V e ρ e φ e i (X ) • φ e j (X ) dΩ (2.31)
the external forces and internal forces vectors components are given by the following relations:

F ext i = ∑ e V e f h d • φ e i (X ) dΩ + ∂ h F V e F h d • φ e i (X ) dS (2.32)
and,

F int i = ∑ e V e
Tr[σ e ε e (N e i (X ))]dΩ (2.33)

In the case of small displacement one can introduce the stiffness matrix K s giving place to the calculation of internal forces through the product of the latter with the displacement vector (2.34)

F int = K s u h (2.34)
with,

K i j = ∑ e V e Tr[K s ε e (N e i (X )) • ε e (N e i (X ))] dΩ (2.35)
In order to establish the numerical governing equations for the solid sub-domain, the physical domain is discretized by using the finite-element method where u s , v s and a s , are the displacement, velocity and acceleration vectors, respectively. In an explicit framework, for the solid sub-domain, matrix M s represents the positive diagonal mass matrix whose components are related to all the DOF of the system. In the following, we assume for simplicity a linear elastic behavior, allowing us to consider a positive definite stiffness matrix K s , which remains constant throughout the calculation. The purpose of the latter assumption is to simplify the first step of the presentation of the solid sub-domain's equations.

Hamilton's principle is used to prescribe equilibrium of the solid sub-domain over a time-interval. Then, the governing equations are integrated in time according to an explicit Central Difference scheme. The current section aims at describing how the coupled solid equations are obtained and how they are made compatible with the formulation given in the previous section.

Following [PRA 04], the equilibrium equations for the coupled solid sub-domain can be obtained by expressing its augmented Lagrangian:

L = 1 2 v T s M s v s - 1 2 u T s K s u s + λ T C s v s + f ext T s u s (2.36)
where f ext s are the external forces to the solid sub-domain. Multipliers λ are used in order to impose the kinematic continuity condition at the interface between the fluid and solid. This condition will be imposed here only onto the normal component of velocity at the interface as it will be subsequently explained.

In (2.36), C s is a constraint matrix that extracts the interface degrees of freedom from the total velocity vector. In general, for matching meshes in geometrically linear structural dynamics, this matrix is binary (containing only zeros and ones), constant and of dimension (N I × N N ), with N I being the number of DOF found at the interface and N N the total number of solid DOF in the sub-domain. However, for fluid-structure sub-domain coupling, we restrict our coupling to normal velocities to the interface. This will reduce the number of DOF along the interface (N I ) by a factor of 2 in 2-D and by a factor of 3 in 3-D. More details about this matrix will be given in the following section. In order to distinguish the case where the constraint matrix is used for FSI coupling along the normal direction only, matrix C s will be referred to as L s . In addition, the interface forces denoted usually by C T s Λ in structural dynamics will be replaced here by L p Λ.

To complete the system of equations, we follow [ROU 90] and [GRA 01] to obtain the definition of the interface energy over the time interval [t 0 ,t e ]:

2. Fluid-solid coupling strategy

L link = t e t 0 λ T N s ∑ k=1 L k v k dt (2.37)
Taking equation (2.37) and the variational of equation (2.36) to apply Hamilton's principle on the interval [t 0 ,t e ] (for further details see [PRA 04]), one obtains in the end:

M s a s + K s u s + L p Λ -f ext s = 0 (2.38
) and:

N s ∑ k L k v k = 0 (2.39)
where we have introduced the Lagrange multipliers Λ. The quantity Λ is homogeneous to a pressure induced by the fluid on the solid sub-domain and it is found by solving the coupled FSI problem.

It is important to note that matrix L p has the same dimension as the matrix L T k but contains information about the surface area of each element along the interface such that when the product L p Λ is performed, the generalized force vector related to the water pressure on the solid is retrieved [LI 15]. It can be seen as an interface force acting on the solid sub-domain.

Furthermore, using the same method, we obtain the velocity constraint at the interface (2.39) expressed for any sub-domain k ∈ [1, N s ], N s being the total number of subdomains. As explained previously, L k are matrices that extract the interface velocity degrees of freedom and project them along the normal direction.

Here we will focus on equation (2.38), which for any time t (with t ∈ [t 0 ,t e ]), guarantees dynamic equilibrium of the solid sub-domain with respect to external forces f ext s and L p Λ.

In case of non-linear behaviour (large displacements and/or non-linear constitutive laws), the internal forces can no longer be represented by K s u s . Instead, the internal forces must be computed by using the following expression:

f int = Ω s B T s σ s dΩ s (2.40)
where Ω s represents the solid domain (which changes in time but is known at the beginning of each time-step since explicit time-integrators are used), B T s is the derivative matrix of the shape-functions in the current configuration and σ s is the Cauchy stress Governing equations for the solid sub-domain tensor.

Taking into account the definitions of f int and f ext s , equation (2.38) can be expressed as:

M s a s + f int = f ext s -L p Λ (2.41)
We will now focus on the time-stepping aspects so as to update the solid configuration from the configuration known at the beginning of the time-step. Let us discretize the time interval [t 0 ,t e ] into equally spaced time-steps ∆t = t n+1 -t n . We can write the equation of motion in (2.41) at t n+1 as:

M s a n+1 s + f int n+1 = f ext n+1 s -L n+1 p Λ n+1 (2.42)
Discretisation of equation (2.42) with the central difference explicit Newmark algorithm gives:

           a n+1 s = [M s ] -1 f ext n+1 s -f int n+1 -L n+1 p Λ n+1 v n+1 s = p v n+1 s + 1 2 ∆t s a n+1 s u n+1 s = p u n+1 s (2.43)
where p v n+1 s and p u n+1 s are called the predictors for the velocity and displacement vectors, respectively, since they depend exclusively on the expression of the kinematic vectors from the known (initial) configuration. The main feature of the explicit timestepping lies in the fact that the new configuration in terms of displacement is known, because it only depends on known quantities at the beginning of the time-step, providing the advantage of being non iterative even in the non linear case contrary to the implicit time-stepping schemes.

The expression of the predictors is given in (2.44).

     p v n+1 s = v n s + 1 2 ∆t s a n s p u n+1 s = u n s + ∆t s v n s + 1 2 ∆t 2 s a n s (2.44)
Finally, equations (2.43) can easily be made compatible with the previously presented fluid equations (2.27) by dividing the time-step by a factor of two in order to obtain the solid configuration at t = n + 1/2. This was previously presented in [ZHE 13].

Coupling strategy

As stated in the previous section, the coupling strategy is based on a dual Schur approach, along the lines of [ZHE 13, GRA 01], by prescribing the kinematic continuity at the interface in terms of normal velocities.

For the present case, we only need to express equation (2.39) for the solid (subscripted s) and the fluid (subscripted f ) as:

L s v s + L f v f = 0 (2.45)
Equation (2.45) prescribes normal velocity equilibrium at the interface between the fluid and the solid.

By guaranteeing that equation (2.45) will be enforced for any t ∈ [t 0 ,t e ], from (2.37), we have:

t e t 0 λ T L s v s + L f v f dt = 0 (2.46)
Considering that the outward-pointing normal vector n s (x k ), for any solid element k, such that x k is positioned at the geometrical barycenter of the interface element (mid-point for a bar, centroid for a triangle), and that it is opposite in sign to the outward pointing normal of the fluid sub-domain (Figure 2.5), the following assertion can be put forward:

n k s (x k ) = -n k f (x k ) = n k (2.47) Coupling strategy
By applying condition (2.47), for an element k, condition (2.45) can be expressed as:

n k • v s (x k ) -v f (x k ) = 0 (2.48)
Expressing equation (2.48) in matrix notation for the whole system and bearing in

mind that v f k (k) = n k • v f (k) for k ∈ (1, N k
), the following system is obtained:

L s v s + v f k = 0 (2.49)
As previously introduced, L s is a geometrical operator enabling the sub-domain DOF to be projected along the interface's normal direction.

For example, in a 2-D configuration, for any interface element k we have:

L s v s (k) = N k 1 v k i j 1 + N k 2 v k i j 2 (2.50)
where the normal vector to element k has been expressed as

N k = [N k 1 , N k 2 ]
assuming k has connectivity to nodes i and j and we have expressed its components in a 2-D Cartesian base (e 1 , e 2 ). We consider from (2.50) that the average velocities of element k in directions 1 or 2 as v k i j 1 or v k i j 2 are defined by:

v k i j 1 = v k i 1 + v k j 1 2 , v k i j 2 = v k i 2 + v k j 2 2 (2.51)
Carrying out the operation described in (2.50) for all surface elements k makes the formulation compatible with the information coming from the fluid side (vector v f k ). Additionally, the dimension of the problem is reduced to (N k × 1).

Remark: From equation (2.38) in the previous section, the product L p Λ allows to link the fluid pressure field to a external nodal force vector f ext f of dimension (N N × 1). The operator L p is very similar to the transpose of operator L s , but includes the value of the surface area s k of element k multiplied to the coefficients of the corresponding column k ∈ (1, N k ). Considering once again a 2-D example, for a single interface element k, after carrying out the operation f ext f = -L p Λ, one gets for each DOF of element k the following expression: One must take now into consideration the time integration aspects of the velocity coupling equation. Since for the fluid integration, we are making use of a two-step Runge-Kutta scheme, we must consider equation (2.49) at times t n+1/2 and t n+1 :

f ext f (I) = - 1 2 N k α s k Λ(k) (2.
   L n+1/2 s v n+1/2 s -v n+1/2 f k = 0 L n+1 s v n+1 s -v n+1 f k = 0 (2.53)
In geometrically non-linear cases, since we can no longer consider the reference geometry for the whole duration of the simulation, operators L s and L p have to be updated at each time-step. For 3-D cases and multi time-step calculations, in which the solid subdomain uses a much smaller time scale than the fluid, the updating of these quantities at each micro time-step can be very taxing on the efficiency of the calculation. One seeks to make the determination of these operators more efficient by choosing to update them at a lower frequency or by storing them using sparse matrix structures and doing operations involving them with dedicated sparse matrix calculation libraries [DEM 99]. In the current work the operators are updated at each time-step since the problems are non-linear.

Contact force modelling

The numerical simulations within the scope of this work will rely on the modelling of contacts between solids in order to predict the damage induced by a flow-driven stone on a functioning hydroelectric turbine blade. Following the explanations put forward in the introductory chapter, where the basic concepts regarding the definition of contacts in a numerical framework and the algorithms for contact detection were discussed, the current section will detail the way in which the contact reaction forces are calculated and how are they made compatible with the fluid-structure coupling effects that were described previously.

Following [BEL 91], we will us the principle of virtual power to express the conservation equation in its variational form. For the sake of simplicity, only two solid bodies will be considered which will be represented by Ω = Ω A ∪ Ω B . The same notations previously used will be used to denote their respective domain for imposed kinematic boundary conditions, ∂ u Ω A and ∂ u Ω B , and the dynamic boundary conditions, ∂ F Ω A and ∂ F Ω B , hence

∂ u Ω = ∂ u Ω A ∪ ∂ u Ω B and ∂ F Ω = ∂ F Ω A ∪ ∂ F Ω B .
Using the principle of virtual power, we can retrieve the variational formulation of the conservation equations for the system comprising both solids :

δP = Ω ∂(δv i ) ∂x j σ ji dΩ - Ω δv i ρb i dΩ - ∂ F Ω δv i f d i dΓ + Ω δv i ρ ∂v i ∂t dΩ = 0 ∀v ∈ U = v(X,t)|v ∈ C 0 (Ω), v = U 0 on ∂ u Ω ∀δv ∈ Û = δv(X)|δv ∈ C 0 (Ω), δv = 0 on ∂ u Ω (2.54)
In equation (2.54), v is the trial velocity vector, δv is a test velocity function that is kinematically admissible on ∂ u Ω, σ is the strain rate tensor, ρ is the solid density and b is the body force vector and f d .

In order to take into account the contact forces, Lagrange multipliers will be used to enforce dynamic continuity along the contact frontier ∂ c Ω between Ω A and Ω B . Following [BEL 91], in order to define the contact problem under a variational form, the trial and test function subspaces in relation with the Lagrange multipliers λ must be defined :

Λ = ∀λ|λ ∈ C -1 (∂ c Ω), λ ≥ 0 Λ = ∀δλ|δλ ∈ C -1 , δλ ≤ 0 (2.55)
where the function in the trial space respect the hypothesis that the forces induced through the contact (Lagrange multipliers λ) can only be of compressive nature (no inter-body adherence). The variational formulation of the initial solid problem plus the variational formulation of the contact problem gives:

δP c = δP + ∂ c Ω δ(λγ n ) dΓ ≥ 0 (2.56)
where γ n was previously defined in Chapter 1 as the normal velocity equilibrium condition along the contact interface.

Discretization of the contact forces

In order to determine the discrete reaction forces induced by the contacts between the solid domains modelled by Lagrangian FEM meshes, this subsection will describe how the contact equation model can be integrated to the semi-discrete governing equation system for the solid described in section 2.2.1.

The use of Lagrange multipliers in order to determine the imposed contact forces on the solid system as well as to calculate the reaction forces of permanent boundary conditions adds an extra unknown to the semi-discrete equations presented in the previous section. This extra unknown calls for the inclusion of an additional algebraic relation between the unknowns of the system. Similarly to the case of the fluid-structure interaction, the additional algebraic relation comes in the form of an added constraint on the velocities. The set of all kinematic constraints (both permanent and contact induced) provide an algebraic system of constraints which for simplicity will be presented as the following equality:

L n+1 c v n+1 s = b n+1 (2.57)
where L c is a geometrical operator allowing to extract the DOF of the system on which contact has been detected through one of the previously presented contact-detecting algorithms and operator b is a vector of known coefficients. Since through the rigid motion of the solid bodies the involved DOF in the contact problem change through time, all of the operators in (2.57) must be recalculated with each new time step.

The Lagrange multiplier representing the reaction forces to the constraints coming from contacts or from imposed boundary conditions appear in the equation of motion (2.41) through an added term:

M s a s + f int = f ext s -L p Λ + L T c Λ c (2.58)
Equation (2.58) is to be solved at each time-step in order to obtain the updated values of the acceleration and thus determine the new solid configuration of the structure on the following time-steps. In order to do so, the internal forces (f int ) and the external forces (f ext s ) need to be obtained as well as the Lagrange multipliers coming from both the constraints (Λ c ) as well as the fluid-structure interaction effects (Λ). In the current study, different specialized solvers will be used to carry out the solid and fluid calculations.

Hence, in an approach mimicking the artificial skin method by Bourel et al. [BOU 06],

the solid solver will retrieve the Lagrange multipliers coming from the contacts separately from the Lagrange multipliers associated to fluid-structure interactions. For this matter, the term related to the fluid-induced forces will be momentarily omitted from the solid equilibrium equation:

M s a s + f int = f ext s + f cons s (2.59)
where f cons s = L T c Λ c represents the forces coming from the contact and permanent constraints prescribed on the solid. Following [CAS 03], equations (2.57) and (2.59) can be used to obtain the algebraic system allowing to retrieve the Lagrange multipliers associated to the constraints on the solid sub-domain:

H n+1 c Λ n+1 c = W n+1 (2.60) with :      H n+1 c = L n+1 c M -1 s L T c n+1 W n+1 = 2 ∆t s (b n+1 -L n+1 c v n+1/2 ) -L n+1 c M -1 s (f ext s n+1 -f int n+1 ) (2.61)
For equation (2.60), operators H n+1 c and W n+1 are calculated from known values at the beginning of the time iteration. The iteration superscript has been explicitly added in order to highlight the fact that when dealing with contacts, these operators have to be calculated at each new iteration (which would not be the case if the only kinematic constraint considered were permanent).

Remark:

The fact that the right hand side of equation 2.60 includes a velocity term calculated at n + 1/2, implies that the total velocity should be known for the concerned degrees of freedom in order to calculate the value of the contact Lagrange multipliers correctly. However, if fluid-structure interaction effects are present on the same degrees of freedom on which the contact is acting, there will be a time-step offset between the effect of the fluid-structure interaction presence and the effect of the contact forces (due to the decoupling-based solving procedure that will be explained in the next chapter). In general, for the current work, there will rarely be cases where a concerned degree of freedom is loaded by contact induced forces as well as by fluid pressure induced forces since, in order for a contact to exist, the fluid separating two potentially contacting solid 2. Fluid-solid coupling strategy surfaces should be evacuated for the contact to be effectively established.

Following the retrieval of the contact Lagrange multipliers Λ c , the contact force vector f cons s can be retrieved and the whole lot of forces acting on the solid system that do not depend on fluid-structure interaction effects can be retrieved. The forces coming from the constraints can be considered to be external factors to the solid sub-domain in question and in turn, they can be regrouped along with the f ext s vector into a more general external force vector f ext tot : As such the more general external force vector can be expressed by :

f ext tot = f ext s + f cons s (2.62)
Using the more general expression the solid equilibrium equation becomes:

M s a s + f int = f ext tot (2.63)
In order to take into account the forces induced by the fluid-structure interaction in the solid, the solid discrete momentum conservation equation becomes:

M s a s + f int = f ext tot + L p Λ (2.64)
Hence, the retrieval of the contact forces is made compatible with the proposed fluidstructure interaction technique that will be described in detail in the following chapter.

Conclusion on FSI couping with the SPH-FEM methods

Over the current chapter, a detailed description of the discretization techniques used to solve and couple the governing equation for the fluid and solid sub-domains has been given. As it was pointed out, the fluid sub-domain will be discretized in space using a variant of the SPH method known as the SPH-ALE method. This technique, in addition to allowing the transport velocity of numerical particles to be different from the velocity of the material continuum, provides the necessary formalism to impose boundary conditions as the solution to a partial Riemann problem. As for time discretization, a explicit Runge-Kutta 2 (midpoint version) scheme is applied. For the solid sub-domain, regarding the spatial discretization, a brief overview of the widely known Finite-Element method has been provided. The time integration scheme used for this sub-domain is the explicit Newmark Central Difference scheme. In order to couple both sub-domains, normal Conclusion on FSI couping with the SPH-FEM methods velocity continuity is imposed at the interface with the use of Lagrange multipliers representing the fluid pressure at the interface. This technique allows for the coupling of different numerical schemes and, as it will be explained in the following chapter, even if different time-steps are used for each sub-domain. Finally, by considering the contact forces between solids to be exclusively a part of the calculations carried out internally by the solid solver (without influence of the coupled calculation), the coupling technique can be made compatible with the contact force modelling approaches present in the solid solver. In the next chapter, the full description of the steps taken to carry out coupled calculations will be elaborated as well as a technique enabling the use of different time-steps in each sub-domain.

Chapter 3 FSI Coupling with different time scales and solid impacts

In this section, the fluid and solid coupling strategy presented in the last section will be extended in order to handle different time scales in each sub-domain.

Incompatible time-step implementation

The use of explicit time integrators on both sub-domains allows for a simplification of the system of equations to solve since the primary variables (displacements for the solid and velocities for the fluid) at the end of the time-step only depend on previously known quantities. This excludes the need for iterative procedures to update the solid and fluid status.

However, explicit integration schemes are known to be prone to instabilities. When using the same time-step in both sub-domains, the efficiency of the coupling algorithm is governed by the size of the smallest time-step in the whole fluid-structure interaction problem. In order to optimize the efficiency and stability of the algorithm, it becomes important to be able to integrate each sub-domain with its own time-step, while respecting the stability of the coupling scheme.

For this study, the use of a smaller time-step is investigated in both the solid and fluid sub-domains. However, the implementation differs for each case. The details regarding the implementation of each case will be elaborated in the current section followed by validation examples given in the following section.

When working with explicit time-integrators, in order to ensure the stability of the calculation, The CFL (Courant-Friedrich-Lewy) condition has to be considered for each sub-domain:

∆t ≤ K CFL × min ∆x i C 0 (3.1)
where ∆t is the required time-step, ∆x i is the characteristic size of the mesh element or SPH particle, respectively, C 0 is the speed of wave propagation in the domain and K CFL ∈ (0, 1). Considering a weakly-compressible fluid state law, for equally sized initial meshes and particle sizes, the application of the CFL condition will require a smaller time-step for the domain where the wave propagation is the fastest. Hence, if the domains considered are water for the fluid and stiff and dense materials for the solid, the smaller time-step is usually related to the solid.

The implementation of smaller time-steps on the fluid side is also an important aspect to consider. Indeed, despite the lower magnitude of the sound speed in the fluid sub-domain, it is important to consider case scenarios where the fluid particle size is much smaller than the element size in the solid, thus producing a much stiffer (in time integration terms) sub-domain. Therefore, work has also been carried out in order to implement a technique Incompatible time-step implementation allowing to integrate the fluid sub-domain with smaller time-steps than in the solid.

3.1.1

The GC method with smaller time-steps in the solid and large interface displacements

Let us consider first the case where the smaller time-step is in the solid sub-domain. Let us now assume that the fluid sub-domain's time-step is a multiple of the solid one:

∆T f = m∆t s , m ∈ N * (3.2)
In this paper, the incompatible GC time-step integration method [GRA 01] is used, because its application to geometrically non-linear interfaces is rather simple to implement. This method enforces velocity continuity at the smallest time scale with the use of Lagrange multipliers.

Remark 1: In the case of explicit time-stepping, the coupling GC method, based on the fine time scale, is more appropriate than coupling techniques based on the coarse time scale because the evolution of the interface geometry and operators is better captured. If the linking of sub-domains were done at the macro time-scale, hypotheses about the Lagrange multipliers and normal vectors to the interface need to be made, leading to a reduced accuracy for the prediction of the interface configuration at the fine time scale.

Remark 2: Nevertheless, when using the GC technique, one must be aware of its inherent dissipative nature which stems out from the interpolation done on the kinematic quantities coming from the sub-domain with the bigger time scale in order to solve the interface problem at the micro time-scale. The fact that the GC method dissipates some energy at the interface when coupling sub-domains with different time-steps can be considered as a drawback for this method. However, three points should be considered: Firstly, it has been observed in many applications [GRA 01, PRA 04, BRU 15] that the method only dissipates a small amount of energy at the interface. Moreover the amount of numerical energy that is lost at the interface can be computed and compared to the overall energy balance as it is usually accepted for Hourglass energy [BEL 84], for example. Secondly, the dissipative features of the GC algorithm have been highlighted when coupling non-dissipative time integrators; the numerical dissipation has been also observed when coupling dissipative time integrators but the amount of dissipation may be much less than the one coming from the time integrator itself or from viscous terms introduced into the formulations. Often, the SPH formulations add viscous terms to stabilize the spa- In this section, the main theoretical aspects behind the implementation of the multistep GC method are reviewed. An algorithm illustrating how the technique is used for fluid-structure interaction simulations will follow. Now let's assume that the large time-step is discretized in m equal sub-steps ∆t s . A typical micro sub-step is denoted j ∈ [1, m]. For the GC method, exact equilibrium is prescribed through the interface at time-step ∆T f = t s 0 ,t s m . Continuity of kinematic quantities is prescribed at any intermediate time j∆t s through an interpolation.

For the macro time-step ∆T f = m∆t s , from (2.24), we have for the fluid sub-domain:

M m f v m f k + Λ m = G m f (3.3)
For the solid, we have for each micro time-step j∆t s :

M s a j s + f int j s = f ext j s -L j p Λ j (3.4)
By combining equation (2.43) with equation (2.49), we can express the velocity constraint condition for any micro-time-step t j as:

L j s p v j s + ∆t s 2 a j s + v j f k = 0 (3.5)
Equations (3.4) and (3.5) can be expressed in matrix form as:

∆t s 2 M s L j p L j s 0 a j s Λ j = G j s W j s (3.6)
where:

G j s = ∆t s 2 (f ext j s -f int j s ) (3.7)
W j represents the constraint condition for any time j∆t s . As in the case featuring the mono time-step coupling procedure [ZHE 13], this condition only involves the normal component of the velocity. For the current case we have:

W j s = -L j s p v j s -v j f k (3.8)
where v j f k is the linear interpolation of the boundary velocity from the fluid side known at j = 0 and j = m:

v j f k = 1 - j m v 0 f k + j m v m f k (3.9)
For the last micro time-step ∆t s (when j = m), we can write the same system of equations:

∆t s 2 M s L m p L m s 0 a m s Λ m = G m s W m s (3.10)
where W m s is no longer an interpolated quantity. Thus exact equilibrium can be prescribed at the macro time-scale.

Regarding the equations from the fluid side, there is no change to what was previously presented since the only unknown needed to update the fluid state is the boundary velocity in the normal direction, which is obtained at the macro-time scale, i.e:

M m f I I 0 v m f k Λ m = G m f W m f (3.11) with W m f = -L m s v m s .
A bordered approach is used to solve this system of equations. This means that the solution to the coupled problem is the superposition of the solution to an unconstrained or free problem and to a constrained or linked problem.

For t = t m one has:

v m s = v f ree m s + v link m s v m f k = v f ree m f k + v link m f k (3.12)
From (3.6) we can deduce the expressions for the free accelerations in the solid at any time t = t j :

a f ree j s = [M s ] -1 G j s (3.13)
Using (3.13), the expression of the free solid velocity is:

v f ree j s = p v 0 s + ∆t s 2 a f ree j s (3.14)
From system (3.11), the expressions of the boundary free velocities of the fluid are 3. FSI Coupling with different time scales and solid impacts directly obtained at t = t 0 and t = t m :

     v f ree 0 f k = M 0 f -1 G 0 f v f ree m f k = M m f -1 G m f (3.15)
Finally, from (3.6) and (3.11), the following expressions concerning the linked velocities in each sub-domain can be deduced:

             v link 0 f k = M 0 f -1 Λ 0 v link m f k = M m f -1 Λ m v link j s = ∆t s 2 [M s ] -1 L j p Λ j = ∆t s 2 a link j s (3.16)
To obtain the constrained and unconstrained boundary velocities from the fluid at any t j , we use a linear interpolation between the values at t 0 and t m :

       v f ree j f k = 1 - j m v f ree 0 f k + j m v f ree m f k v link j f k = 1 - j m v link 0 f k + j m v link m f k (3.17)
We can now express equation (2.49) at any t j in terms of its constrained and unconstrained values:

v link j f k + L j s v link j s = -v f ree j f k -L j s v f ree j s (3.18)
where the right-hand-side of equation ( 3.18) will be referred to as -b j . Notice that vector b j exclusively depends on the unconstrained quantities coming from the fluid and solid sub-domains. Hence, we express (3.18) as:

v link j f k + L j s v link j s = -b j (3.19)
Inserting the expressions (3.16) and (3.17) into (3.19), the following expression is obtained:

- 1 - j m M 0 f -1 Λ 0 + j m M m f -1 Λ m - ∆t s 2 L j s [M s ] -1 L j p Λ j = -b j (3.20)
In order to apply the GC method, we approach our coupling formulation to what was established in [GRA 01] for the coupling of solid sub-domains. To do this, we begin by considering a constant M f matrix over the macro time-step:

M f ≡ M 0 f ≡ M m f (3.21)
In practice, this matrix can be taken equal to its value at t = t 0 or t = t m or to an average value between both. Here, we have decided to use a constant value for the M f matrix over the macro time-step, which we have taken equal to M m f . Additionally, a hypothesis regarding the value of the Lagrange multipliers at the micro time-steps is made by considering once again a linear interpolation between the values at at t = t 0 and at t = t m :

Λ j = 1 - j m Λ 0 + j m Λ m (3.22)
By inserting equations (3.21) and (3.22) and factoring out the Lagrange multiplier expressed at t j , the following relationship is established:

M f -1 + ∆t s 2 L j s [M s ] -1 L j p Λ j = b j (3.23)
The operator inside the brackets pre-multiplying Λ j is called the condensed operator H j .

The system becomes:

H j Λ j = b j (3.24)
and can be solved for at each micro-iteration j.

Remark: By retrieving the Lagrange multiplier, the value of the link velocities can be obtained by using (3.16) and thus the updated configuration for the solid U j s for any t = t j . If t = t m the value of the boundary velocity is sent back from the coupler to the fluid solver and the new fluid configuration can be obtained (U m f → U 2m f ). Figure 3.1 schematically presents the coupling procedure implemented for the integration of both the fluid and the solid sub-domains from t 0 to t m . The multi-step coupling strategy can easily be adapted to take into account the mid-point stage needed by the fluid's time integration. The macro time-step is divided in two and the procedure is done twice. Finally, the procedure used over the macro time-step for the multi-step coupling algorithm is detailed in Algorithm 1. The flowchart in Figure 3.2 describes the way in which the algorithm is implemented. Colors and acronyms have been adopted to describe the tasks done by the different solvers involved in the coupling procedure: SLS represents the solid solver Europlexus , FLS stands for the fluid solver ASPHODEL and CPL is the internally developed coupling software.

s s s s s s s C C C C C C C F F F t j =t 0 t j+1 t j+2 t j =t m
In practice, the position of the solid interface at time t j = t m = t n+1/2 or t j = t m = t n is not known exactly as it would be the case for m = 1 (same time-step in both domains). Once the solid interface (position and normals) is known, the fluid solver determines the interface elements under FSI. One must therefore first provide an approximative interface status Ĩs to the fluid solver. For simplicity in the implementation of the current algorithm, it has been decided to use the known position of the wall at the first micro time-step j = 1 as the estimated interface status Ĩs . This position is then sent to the fluid solver in order to get back a list of wet surface solid elements. This choice was adopted because the list of elements under FSI is nearly constant between two macro time-steps for the cases performed up to now as the solid motion is rather slow. Hence, the interface geometry change is very small. Moreover the fluid time-step is rather small and the fluid position change is very limited within each macro time-step. Additionally, the time integration within the fluid sub-domain is only performed at j = m, hence, the precise position of solid interface at intermediate times is not used.

for n = 1, n ≤ nMax do -FLS obtains U m f from U 0 f and v 0 f k for j = 1, j ≤ m do SLS computes u j s = p u j s → f int j (u j s ), f ext j if j = 1 then
-SLS calculates predicted interface status Ĩs and sends it to FLS → FSI exists?

if FSI = true then -FLS computes M m f and G m f and creates ID list of FSI elements, m ID -FLS sends M m f , G m f to CPL and m ID to SLS, waits until j = m -SLS calculates a f ree 1 s , v f ree 1 s , L s 1 and L p 1 -SLS sends L s 1 v f ree 1 s and L s 1 M s L p 1 to CPL -CPL estimates interpolations of M 1 f and G 1 f and uses them to estimate v f ree 1 f k -CPL calculates H 1 and b 1 to obtain Λ 1 and sends it to SLS -SLS calculates a link 1 s and v link 1 s → obtains a 1 s and v 1 s → new configuration U 1 s else -SLS obtains a 1 s and v 1 s → new configuration U 1 s end if end if if j > 1 and j < m then if FSI = true then -SLS calculates a f ree j s , v f ree j s
, L s j and L p j -SLS sends L s j v f ree j s and L s j M s L p j to CPL

-CPL calculates v f ree j f k through M 0 f , G 0 f , M m f , G m f -CPL calculates
end if end if if j = m then if FSI = true then -SLS calculates a f reem s , v f reem s , L sm and L pm -SLS sends L sm v f reem s and L sm M s L pm to CPL -CPL calculates H m and b m to obtain Λ m and sends it to SLS -CPL calculates v linkm f k , hence v m f k -SLS calculates a linkm s and v linkm s → obtains a m s , v m s → new configuration U m s -SLS

Incompatible time-step implementation

Notice that for Algorithm 1 is given for a general time macro time-step ∆T f , regardless of the fluid sub-stage considered, i.e. t m = t n+1/2 or t m = t n+1 . One should nevertheless recall that the solid kinematic quantities are not updated for the first stage of the calculation (t m = t n+1/2 ) as this stage only provides interface information that the fluid needs in order to advance in time.

Remark: Several possibilities can be considered to estimate the interface position: One could calculate an estimation of the position of the wall using only uncoupled kinematic quantities and those known at j = 0. Alternatively, one could iterate with the fluid solver in order to update the list of FSI elements at each micro time-step as the exact position of the wall is known for the next micro time-step at the start of each micro time-step: This last method is rather heavy and CPU time consuming.

Smaller time-steps in the fluid

Despite wave propagation being usually faster within the solid, cases concerning highly dense liquids or of smaller particle size should be considered. These cases call for the implementation of a smaller time-step within the fluid. The GC technique can be applied as well, however the implementation slightly differs. Since the solid sub-domain is being integrated with a larger time-step, the position of the wall at any instant t j can be directly interpolated between the known positions at t 0 and t m :

u j f k = (1 - j m )u 0 f k + j m u m f k (3.25)
Knowing the solid wall's position at each micro time-step, the fluid solver can determine the FSI element list at any time t j under consideration. The information regarding the elements under FSI as well as the coefficients M j f and G j f needed for the coupling are then transferred to the coupling software in order to determine the updated coupling constraint matrices L j s and L j p . For this coupling procedure, more information is transferred to the coupler in order to reduce the tasks done by the solid solver and improve the efficiency of computation of the coupled system.

With the implementation of Algorithm 2, the velocity at the interface for any t j can be obtained:

v j f k = v f ree j f k + v link j f k (3.26)

FSI Coupling with different time scales and solid impacts

Knowing the interface velocity at time t j (3.26), the fluid solver can complete its time integration up to the next time-step, thus obtain the new fluid configuration U j+1 f . The Runge-Kutta 2 mid-point time-step is still used in the fluid to go from time t = t j to t = t j+1 . Since the explicit Newmark scheme used for the solid is second order accurate, there is no need to carry out a solid calculation at the mid-step.

Finally when t = t m the coupling software sends back the result of the coupled calculation (the Lagrange multiplier vector), for the solid to complete its integration at the new configuration U m s . Figure 3.3 schematically shows how the integration procedure is carried out. Some tasks that were done previously by the solid solver when using smaller timesteps in the solid, are now carried out by the coupling software in order to avoid the need to make use of the solid solver at instants that are intermediate to the coarse time-scale. For this purpose, the solid mass matrix M s , which is diagonal and invariant, as well as the connectivity information regarding the boundary elements is transferred at the beginning of the calculation to the coupling software. Additionally, at the beginning of the large time-step related to the solid, the position of the wall at times t 0 and t m as well as the free velocities of the solid at t 0 and t m , v f ree 0 s and v f ree m s

s s s C C C C C C C F F F F F F F t j =t 0 t j+1/2 t j+1 t j =t m
, respectively, are transferred from the solid solver to the coupling software.

The value of the unconstrained or free velocity from the solid side can be interpolated to obtain its value at t = t j :

v f ree j s = (1 - j m )v f ree 0 s + j m v f ree m s (3.27)
From equation (3.11), the free velocity at the fluid interface is:

Incompatible time-step implementation v f ree j f k = M j f -1 G j f (3.28)
Using (3.28), the right hand-side vector b j of the interface equation can now be obtained:

b j = v f ree j f k + L s j v f ree j s (3.29)
The interface equation can thus be written as:

-M j f -1 Λ j -L j s v link j s = -b j (3.30)
By adapting equation (3.16) to the current case and by using a linear interpolation for the solid link velocities at times t = t 0 and t = t m to find its value at t = t j , equation (3.30) can be expressed as:

-M j f -1 Λ j -L j s (1 - j m ) ∆t s 2 [M s ] -1 L 0 p Λ 0 + ( j m ) ∆t s 2 [M s ] -1 L m p Λ m = -b j (3.31)
where the operators L j s and L j p are obtained by the coupling software. Finally in order to establish the interface equation for the case where the smaller timesteps are found in the fluid sub-domain, the following hypothesis is considered

L j p Λ j = (1 - j m )L 0 p Λ 0 + ( j m )L m p Λ m (3.32)
By using (3.32), (3.31) can be expressed as (3.24), with:

H j = M j f -1 + ∆t s 2 L j s [M s ] -1 L j p (3.33)
Finally, by solving (3.24), the Lagrange multipliers Λ j can be retrieved and the linked fluid interface velocity at t j :

v link j f k = M j f -1 Λ j (3.34)
The flowchart given in Figure 3.4 gives an overview of how the algorithm is implemented for the case when the fluid sub-domain is integrated using smaller time-steps. For further details, Algorithm 2 is given to detail the implementation of the coupling algorithm.

FSI Coupling with different time scales and solid impacts

Remark: The approach described in this section has only been tested in 1-D and provides accurate results when compared to results obtained when both sub-domains were integrated with the same time-step. However, since the interface problem size is reduced to a single DOF, the volume of information transferred from the solid solver to the coupler and fluid solver is very small. Further testing in 2-D and 3-D should be carried out in order to explore more deeply the quality of the described procedure by studying how the additional data exchanges hinder the efficiency of the method.

for n = 1, n ≤ nMax do for j = 1, j ≤ m do if j = 1 and stage = 1 then -SLS sends [u 0 f k , u m f k , v f ree 0 sk , v f reem sk ] to CPL -CPL sends [v 0 f k , u 0 f k , u m f k ] to FLS -FLS determines if FSI exists if FSI = true then -FLS uses u 0 f k and v 0 f k to obtain U j+ 1 2 f else -FLS directly obtains U j+ 1 2 f end if end if if [ j > 1 and j < m] or [ j = 1 and stage > 1] then if stage = 1 then -FLS obtains u j f k through interpolation and determines if FSI exists if FSI = true then -FLS sends M j f , G j f and list of FSI elements (m ID ) to CPL -CPL interpolates u m f k ,v f reem sk → u j f k ,v f ree j sk → estimates L s j , L p j -CPL calculates v f ree j f k , b j and H j -CPL solves interface problem to obtain j → v link j f k → v j f k -CPL sends v j f k to FLS → U j+ 1 2 f else -FLS directly obtains U j+ 1 2 f end if end if if stage = 2 then -FLS obtains u j+ 1 2 f k → FSI exists? if FSI = true then -FLS sends M j+ 1 2 f and G j+ 1 2 f to CPL -CPL obtains u j+ 1 2 f k ,v f ree j + 1 2 sk → L s j+ 1 2 , L p j+ 1 2 -CPL calculates v f ree j + 1 2 f k , b j+ 1 2 and H j+ 1 2 -CPL solves interface problem to obtain j+ 1 2 → v link j f k → v j+ 1 2 f k -CPL sends v j+ 1 2 f k to FLS → U j+1 f else -FLS directly obtains U j+1 f end if end if end if end for -FLS sends M m f , G m f and m ID to CPL if FSI = true then

Stability of the algorithm

The current section will focus on studying the stability of the coupling algorithm when applied to the case when the smaller time-step is used in the solid. Since the GC technique is also used for the case when the smaller time-step is used in the fluid, a similar stability study can be done for that case that will not be carried out here for the sake of brevity. This study will be carried out by using the energy method [HUG 79] along with the following notations:

x 0 = x m -x 0 (3.35)
and:

x j = x j+1 -x j (3.36)
where x can be any tensor quantity named in the present article. The total energy expression associated with each sub-domain is given by the following expression:

T i (a 0 ) + V i (v 0 ) = -D i ( a 0 ) -E i link ( v 0 , Λ 0 ) (3.37)
where i represents the different sub-domains under consideration with T i being the term associated with the kinetic energy, V i the term associated with the internal energy, D i the one associated with the numerical damping and E i link , the term associated with the interface forces which comes from the coupling algorithm. In order to produce a stable algorithm, one must ensure that this last term is equal to zero or at least smaller than zero. In the latter case, stability is ensured but numerical dissipation occurs, altering the accuracy features of the coupling algorithm with respect to the different time integrators involved in the co-simulation.

Since for the solid we use an explicit central difference scheme with the FEM technique, we can infer that the time and spatial integration of the solid sub-domain neither produces nor dissipates numerical energy [COM 02].

On the fluid side, we use a Runge-Kutta 2 explicit time scheme which is widely known to be stable provided that the time-step size conditions are met. However, as stated previously, numerical dissipation can occur through the spatial integration scheme used in this domain.

Following [GRA 01], the last term in equation (3.37) can be expressed as:

E i link ( v i 0 , Λ i 0 ) = 1 m∆t s v 0 f k T Λ 0 + m ∑ j=1 1 ∆t s L j-1 s v j-1 s T Λ j-1 (3.38)
Using equation (3.12) to split the velocity values into their free and link components, equation (3.38) can be rewritten as:

E i link ( v i 0 , Λ i 0 ) = 1 m∆t s v f ree 0 f k T Λ 0 + 1 m∆t s v link 0 f k T Λ 0 + m ∑ j=1 1 ∆t s L j-1 s v f ree j-1 s T + L j-1 s v link j-1 s T Λ j-1 (3.39)
By using:

Λ 0 = m ∑ j=1 Λ j-1 (3.40)
and the following equation (obtained by expressing v f ree j f k as an interpolation of its values at t = t 0 and t = t m ):

v f ree j-1 f k = 1 m v f ree 0 f k (3.41)
equation (3.39) can be expressed as:

E i link ( v i 0 , Λ i 0 ) = 1 m∆t s v link 0 f k T Λ 0 + m ∑ j=1 1 ∆t s L j-1 s v link j-1 s T Λ j-1 + m ∑ j=1 1 ∆t s v f ree j-1 f k T + L j-1 s v f ree j-1 s T Λ j-1 (3.
E i link ( v i 0 , Λ i 0 ) = 1 m∆t s v link 0 f k T Λ 0 + m ∑ j=1 1 ∆t s L j-1 s v link j-1 s T Λ j-1 - m ∑ j=1 1 ∆t s v link j-1 f k T + L j-1 s v link j-1 s T Λ j-1 (3.43)
where two terms cancel out to give in the end :

E i link ( v i 0 , Λ i 0 ) = 1 m∆t s v link 0 f k T Λ 0 - m ∑ j=1 1 ∆t s v link j-1 f k T Λ j-1 (3.44)
Using equation (3.16), one can express (3.44) as:

E i link ( v i 0 , Λ i 0 ) = 1 m∆t s v link 0 f k T M 0 f v link 0 f k - m ∑ j=1 1 ∆t s v link j-1 f k T M j-1 f v link j-1 f k (3.45)
With (3.21) from section 3.1.1, equation (3.45) can be expressed as:

E i link ( v i 0 , Λ i 0 ) = 1 m∆t s v link 0 f k T M f v link 0 f k - m ∑ j=1 1 ∆t s v link j-1 f k T M f v link j-1 f k (3.46)
Following the lines of [GRA 01] and [COM 02], it can be shown that expression (3.46) can be written as a sum of negative squares by considering that M f is positive and definite. These implications lead to conclude that the multi time-step coupling algorithm induces dissipation at the interface, hence guaranteeing its stability.

Ẽi link = (α n+m -α n ) 2 + m m ∑ j=1 (α n+ j -α n+ j-1 ) 2 (3.47)
One can show that such an expression expands to a sum of negative squares [COM 02].

Through (2.25), M f is observed to be positive and definite and thus the interface velocity vectors are bounded quantities [BEL 73]. These implications lead to conclude that the multi-time step coupling algorithm induces energy dissipation at the interface, hence guaranteeing its stability.

Numerical dissipation at the interface depends mostly on the shape of the velocity curve that is being linearly interpolated. If the shape of the curve differs largely from a linear evolution, the quality of interpolations done will become worst and hence more energy will be lost. There have been several studies done in the past [GRA 01, BRU 15] where the amount of energy dissipated has been showed to be minimal compared to energy dissipated by the employed time integration and space integration schemes.

FSI Coupling with different time scales and solid impacts

Conclusion

A heterogeneous time-step coupling technique has been proposed for fluid structure interaction problems, allowing to carry out accurate and stable FSI simulations. If explicit time integrators in both physical domains are used, the involved system of equations become much simpler to solve. However, the dependence on the time-step size becomes a major drawback when a solid fine mesh or a small fluid particle size were adopted. In order to prevent the time-step size requirements of one sub-domain from being inherited by the other sub-domain, one must be able to integrate each domain with a different timestep. This objective was accomplished by adapting techniques coming from the coupling of solid sub-domains in structural dynamics to fluid-structure interaction problems. Numerical case-studies will be investigated in the following chapter in order to highlight the validity and relevance of the implemented techniques.

Numerical examples

In this section, several numerical validation examples are presented most of which come from [NUN 17].

The first examples will be simple validation FSI cases putting forward comparisons between the multi-step coupling strategy and the previously developed coupling strategy presented in [ZHE 13, LI 15]. The advantages and drawbacks of the multi-step coupling strategy will be highlighted notably by showing the time gain and accuracy of using the proposed strategy but also by showing convergence rate issues and slight dissipation at the interface stemming from the hypotheses made to put the current strategy into use. Following this initial validation examples, cases taking into account the impacts of projectiles inside jet flows will be put forward in order to numerically characterize the damage induced on turbine blades through these type of phenomena. The initial length of the beam is L 0 s = 1 m, its initial solid density ρ 0 s = 2700 kg/m 3 , and its initial section area A 0 s = 0.01 m 2 . Young's modulus is E s = 67.5 GPa. The solid bar is discretized with 100 linear truss finite-elements. The tube of water has also a length of L 0 f = 1 m but contains ten times more particles than elements in the bar. A uniform pressure step of 20 MPa is imposed at the time t = 0 s in the fluid cavity. We will make use of the GC method to integrate the solid sub-domain and fluid sub-domains with different time-steps. The results will be compared with those obtained by [ZHE 13] when using the same time-step to carry out the coupling calculations. 2) shows the comparison between the results of the simulation using the method described in the current paper and another one in which both sub-domains are integrated with the same time-step, i.e. ∆t f = ∆t s = 10 -9 s .

The results show good agreement when using a much finer discretization in the fluid sub-domain leading to a much smaller time-step in the fluid than in the solid in order to satisfy the CFL condition.

Next, for the same test case, we will make use of the GC method to integrate the solid sub-domain with a smaller time-step than the fluid. To perform this simulation 500 truss elements are used to discretize the bar and 100 fluid particles for the tube. The ratio between time-steps, m = ∆t f /∆t s , is equal to 100 with ∆t s = 10 -8 . Once again, we see that good agreement also exists when using a smaller time-step to integrate the solid through the use of the GC method for time-steps that differ considerably.

However, it has to be noted that the interface remains quite small when compared to the size of the whole problem (1 DOF for the interface compared to 500 solid elements + 100 fluid particles).

In order to verify that coupling on the normal component of the velocities is done at

Accuracy study of the proposed algorithm

For the current test-case, we will carry out a time convergence study of the multi timestep coupling algorithm. The strategy proposed in [MIC 04] to study convergence in time provides a simple yet effective way to carry out this study and we will apply it here. The said strategy consists in choosing a starting time-step size τ and carrying out a simulation. Then, simulations using time-step sizes corresponding to 2τ and 4τ are carried out subsequently. To analyse convergence in time, we track the value of a field on a point in the structure or in the fluid (referred to as X num ) and calculate the relative error using an L n norm (with n = 2 in this case) between the simulations using the different time-step sizes. To obtain the order of accuracy in time, denoted by Ω, we compare the natural logarithm of the ratio of subsequent calculation errors to the natural logarithm of 2, i.e.:

Ω = ln X 4τ num -X 2τ num 2 X 2τ num -X τ num 2 ln(2) (4.1)
For this study, we track the displacement of the end of the bar that is in contact with the fluid (u s (L 0 s )). Two studies will be done using different reference time-steps τ, a coarser one at τ f 1 = 1.0 × 10 -6 s and a finer one at τ f 2 = 2.5 × 10 -7 s (for the fluid sub-domain). A smaller time-step will be used in the solid with a ratio of m = 100 between both subdomains. Thus, for the solid we have τ s 1 = 1.0 × 10 -8 s and τ s 2 = 2.5 × 10 -9 s. Table 4.1 gives the results of the convergence study. As it can be seen from Table 4.1, the order of accuracy of the multi time-step coupling algorithm is not second order as it would have been the case when using a mono timestep coupling strategy. As explained in section 3.2, the interpolations and hypotheses done when using the multi time-step coupling strategy induce numerical energy dissipation at the interface, hence the global order of convergence of the algorithm is degraded.

As for the case when the smaller time-steps are used in the fluid, numerical tests carried out for τ f 2 = 2.5 × 10 -7 s using a time-step ratio of m = 10 showed once again a drop in the convergence rate of the algorithm. The results of these latter numerical tests are presented in Table 4.2.

m = 1 Ω(u s (L 0 s )) 2.044 m = 10 Ω(u s (L 0 s )) 1.025
For the following applications, we have considered only the case where the smaller time-step is used for the solid sub-domain as this case is the main target motivating this study.

2-D hydrostatic water on a linear elastic plate

Next, we consider a simple 2-D test case in which we couple a linear elastic plate with a column of water which is initially in hydrostatic equilibrium. Figure 4.5 shows the configuration of this test case: a rigid water reservoir has a geometrically linear elastic bottom which is clamped at both ends, a pressure gradient within the water reservoir is present due to the gravity effect and the fluid pressure at the free surface is taken equal to zero. The other simulation parameters are given in Table 4. 3. particle and mesh size. The speed of sound being roughly three to four times larger in the structure than in the fluid, we will nevertheless use a time-step of ∆t f = 5 × 10 -7 s to integrate the fluid whereas the time-step in the solid will be ten times smaller (m = 10).

The results of the simulation are given in Figure 4.6 in which the displacement of a node located at the middle of the plate predicted by the proposed multi time-step approach is compared to the reference results obtained from an homogeneous time-step approach.

Once again, it can be seen that a good agreement is achieved between the two results.

When running simulations for multidimensional cases, integrating each sub-domain with different time-steps, while ensuring a good accuracy, is crucial for reducing the computation times. Indeed, the fluid calculation is much more time-consuming than the solid calculations. So being able to adopt a coarse time-step for the fluid computations leads to an improved performance.

Breaking dam flow on an elastic wall

Next, another 2-D test case of a fluid-structure interaction problem is considered, whose initial configuration is shown in Figure 4.7. In a rigid wall container, a column of water is initially located at the left side of a container. The water column is in hydrostatic equilibrium. An elastic wall is clamped placed to the right at the middle of the container.

The geometric and discretization parameters are given in Table 4. 4.

The material properties of the solid are such that the initial solid density is ρ 0 s = Once again, due to the same spatial discretization used on both sub-domains, we use a time-step of ∆t f = 2 × 10 -5 s for the fluid whereas a time-step that is 10 times smaller veloped by EDF (Electricité de France).

As presented in [LI 15], there is still a discrepancy between the results obtained in [WAL 05] and those obtained through the use of the present method. However, the results obtained by [IDE 08] come quite close to what was obtained with the proposed method, especially for the coupling done with the Europlexus software. Good agreement between the results obtained when integrating with same and different time-steps is also found, despite the large displacements undergone by the structure. Thus, the interest of the multi time-step GC algorithm for this case study is highlighted. For the last 2-D example, a water jet impacting a steel plate will be simulated. In the current simulation, the water jet will impact the stationary plate with a downward velocity of v f = -100 m/s. Three subsequent calculations will be carried out on which the solid and fluid mesh density will evolve in order to asses the proposed strategy's convergence in time and space. Figure 4.10 presents the geometrical configuration of the numerical 4. Numerical Simulations simulation currently being considered.

L f h f L s h s v f Ω f Ω s M Figure 4
.10: Configuration for the 2-D jet impact on plate case study Concerning the fluid parameters, liquid water values at room temperature will be used for the simulation, hence the reference density ρ 0 f is 1000 kg/m 3 and the reference speed of sound c 0 f is 1500 m/s. The solid plate is considered to be made out of steel, hence the typical material parameters for an isotropic elastic steel will be used. As such the solid density ρ s is 7800 kg/m 3 , the Young's modulus E s is 200 GPa and Poisson's ratio ν is 0.34. Table 4.5 contains the value of the geometrical parameters for the present study: The simulations will be carried out three times with three different sizes for the mesh and particle discretizations. The time-step considered for the solid sub-domain will be the one corresponding to the smallest spatial discretization obtained through the CFL condition. The fluid time step will be 5 times larger than the solid one. The subsequent mesh and particle sizes along with the time-step use for each sub-domain are given in Table 4.6. follows the same tendency as the mesh density is refined. The mesh convergence of the multi-step coupling strategy is thus well represented in the current study and serves as a preamble to the fluid-structure and impact simulations that will be presented next.

Numerical Simulations

FSI simulations with impacts

In the current section, impacts of flow driven projectiles on stationary targets will be carried out. Initially a 2-D parametrical study will aim to show the influence of the size of the projectile on the potential damage that the stationary surface may go through after the impact. Next, 3-D simulations using a stationary plate and followed by a real turbine blade model as target will be studied.

2-D impact simulations, influence of the size of the projectile

For the first study involving impacts and fluid-structure interaction effects, we will consider the case of a circular projectile impacting a flat 2-D plate that is simply supported on its sides and bottom. The projectile engulfed in the water jet will hit the immobile plate with a velocity of v f = -100 m/s. In fact the configuration of the current test case is rather similar to the configuration of the case described in 4.1.4 but with the addition of the contact problem and different material properties and geometrical parameters.

Concerning the simulation parameters for the projectile, some hypotheses will be made in order to create a realistic but computationally efficient scenario. Firstly the material properties of the projectile will be considered to be homogeneous throughout its whole extent.

Since the composition of rocks and gravels found in streams can be quite heterogeneous event within a single specimen, this hypothesis is rather strong but allows to lighten up the computational cost of the whole simulation. Furthermore, in order to position the study in a unfavourable situation for the target body, the projectile will be considered made of quartz which is one of the hardest materials that can be found in river and streams usually mixed into granite composites.

To model the quartz projectile, a Drucker-Prager [DRU 52] elasto-plastic material model will be used which is a commonly used to model rock, concrete and other pressuredependent materials. This material model allows to establish different compression and traction limit stresses after which the material reaches a plastic state. In the solid solver used for the current work, one can even establish a stress limit triggering brittle fracture once a certain stress state is surpassed. Finally in order to simplify simulation time and avoid complicated mesh structures, the geometry of the projectile will be simply circular.

The steel target plate will be modelled using a von Mises elasto-plastic material model which is commonplace for the simulation of metals and alloys. In order to simulate erosion, a post-calculation criterion based on the total plastic strain will be used to determ- The material properties for the rock projectile are its density ρ r = 2200 kg/m 3 , its Young's modulus E r = 75 GPa, its Poisson's coefficient ν = 0.3, its traction stress limit σ T = 50MPa and its compression stress limit σ C = 1.15 GPa.

For the steel plate, the basic material properties include its density ρ p = 7800 kg/m 3 , its Young's modulus E p = 200 GPa, Poisson's coefficient ν = 0.33 and the elastic limit σ e = 200 MPa. After the elastic limit is reached, the non-linear relation between the stress and the strain in the material follows a curve obtained experimentally for this material. A maximum total plastic equivalent strain can be used to determine if erosion is present within the impact zone. For this matter, elements surpassing 8% of total plastic will be considered as eroded [KIM 15, SEI 16]. .13 presents a comparative result of the damage undergone by the steel plate both in the case when the rock impact occurs in a "dry" configuration (top row) as well as in a fully coupled configuration (bottom row). The results presented in this Figure correspond to the case when the radius of the rock is r s = 1 cm. Quantitatively, the damage can be analysed by looking at the total equivalent plastic strain variable given by the solid solver for the target structure. Figure 4.13 exhibits how the damage undergone by the steel plate in the dry configuration is significantly more pronounced than the damage undergone by the plate in the fully coupled configuration. As it can be seen in the first image of the coupled calculation on Figure 4.13, much of the stone's kinetic energy is absorbed by the water jet as the fluid gets trapped between the rock and the steel plate.

Hence, the impact velocity of the rock becomes much lower in the case of the coupled calculation compared to the dry impact case. For the case of the coupled impact, the velocity of impact has been recorded to be around v r ≈ -40 m/s as opposed to the dry case were the velocity of the projectile is the same as the initial velocity of the water jet. In order to take into account the influence of the size of the stone on the extent of the impact damage, a parametric study using three different rock sizes has been carried out where the radii of the circular stones will take the values given in Table 4.7. Figure 4.14 presents a study of the influence of the size of the rock projectile on the damage induced on the target steel plate. The top row corresponds to the case where the projectile has a radius r s = 5 mm, the middle row corresponds to the case where r s = 7.5 mm and finally the bottom row is for the case where r s = 1 cm. For the current study, once again, the total equivalent plastic strain resulting after the impact is taken as the measure of the extent of the damage the plate has gone through. As it is clear from Figure 4.14, the damage induced on the plate is proportional to the size of the projectile. The largest projectile has naturally the most momentum and the water "cushion" forming before impact between the projectile and the plate is less successful at slowing the rock before impact. For the smallest rock, the impact velocity is v s ≈ -19 m/s, for the medium sized one the impact velocity is v s ≈ -30 m/s and finally, as already stated, for the largest rock, the impact velocity is v s ≈ -40 m/s.

Multi-time step coupling strategy gains

The simulations done in the current study were done with a time-step ratio of m = 10 between the fluid and solid sub-domain in order to satisfy the stability time step conditions for each sub-domain. Thus, considering the spatial mesh spacing and the values of the speed of sound inherent to each sub-domain, for the solid we used a time-step of ∆t s = 1.0 × 10 -8 s and for the fluid a time step of ∆t f = 1.0 × 10 -7 s. These simulations were carried out for the case when the radius of the stone is r s = 1 cm. .

To highlight the accuracy of the proposed method even in such a non-linear case as the present one, the vertical displacement of point M, whose position is shown in 4.12, will be tracked in time in and compared to the calculation carried out when both of the sub-domains used the same time-step. The comparative graph showing the results of the mono-time step calculation to those of the multi time-step one are presented in Fig. 4.15. As seen in Figure 4.15, there is good agreement between the calculation done with the same time-step and the one done using a ratio of m = 10 which allowed each subdomain to be integrated with a time step close to its stability time step imposed by the CFL condition. This example highlights the applicability of the proposed method even in the event of large displacements, contacts and material non-linearities.

In order to quantify the time gain when using different time steps for each sub-domain, the total elapsed time of calculations has been recorded for the mono time-step and the multi time-step cases. The calculations presented here were all executed on the same 4. Numerical Simulations server, composed of current generation CPUs clocked at 2.50 GHz. For these calculations, a single CPU was attributed for each solver as well as for the coupling software. Table 4.8 compares the time elapsed for the coupled calculations using the same and different time-steps. As it can be seen from table 4.8, using the multi-time step strategy, allows to reduce the calculation time considerably. For the present 2-D cases, the time-gain might not be a determinant factor for opting to use the multi-step coupling strategy with regards to using the same time step. However, for 3-D cases that will follow, the multi-time step coupling strategy will prove to be crucial to providing results within a reasonable amount of calculation time.

3-D impact simulation

For the first 3-D impact problem presented here, a similar configuration to the one used for the 2-D impact problem will be employed, see Figure 4.12. The stone projectile is naturally replaced by a spheric projectile with radius r s and the rectangular steel target plate is replaced by a rectangular prism. The exact value of the geometrical parameters for the current case study are given in Table 4.9. The material properties of the projectile and the plate will be the same as the ones used in the 2-D case. The same applies for the fluid reference parameters. The speed of the water jet and projectile is once again v j = 100 m/s.

The results of the 3-D impact case are given in Figure 4.16. In Figure 4.16, once again the top row displays the case where the impact is done without coupling and the bottom line the fully-coupled case. As in the 2-D case the damage on the steel plate, quantified by the total elastic plastic strain rate, is higher for the case where no fluid is present. For the coupled case, the fluid absorbs some of the kinetic energy of the projectile with the impact speed being v s ≈ -50m/s with respect to the initial jet and rock velocity of v j = -100 m/s. It is also evident that in the present coupled simulation the damage undergone by the plate is greater than in the 2-D case. The reason for the increased impact velocity in the 3-D case can be explained in two ways.

Firstly the size of the projectile is twice as large than in the 2-D case. The reason for this was mainly to see the effects of a larger projectile on the damage caused by the impact. The second reason for this increased damage can likely come from the fact that in 3-D, the water trapped between the rock and the plate is less constrained than in the 2-D case as it has more direction in space along which it can evacuate. Actually, the simulations done in the 2-D case are not really comparable to the simulations in 3-D presented here since no axisymmetric conditions were used for the 2-D case.

3-D impact simulation on a Pelton bucket

For the last simulation presented in the current work, a more realistic geometric model for the turbine blade as well as for the stone projectile will be used. The material properties of the solid as well as the reference values for the fluid will remain the same as for the case studies presented in subsections 4.1.4 and 4.2.2. The geometrical dimensions of the Pelton turbine bucket are described in Figure 4.17.

In Table 4.10, we can see the detailed dimensions of the Pelton turbine bucket as well as those regarding the water jet and the stone projectile. The particle size for the fluid SPH mesh is given as well and the solid is meshed in a way that the elements that are in contact with the fluid are at least of the same size. The rock projectile has an amorphous shape so no radius can be given as in the precedent case study. Nevertheless, the maximum length in every direction for the projectile is nearly the same and this parameter has been given as the length of the projectile L r .

Using the parameters given, the fluid domain is discritized using 90 000 SPH particles and the solid mesh contains around 86 000 tetrahedron linear elements. The reason why tetrahedron elements are being used is mainly for computational cost reduction. In terms of boundary conditions, the right end of the Pelton bucket in Figure 4.17, is considered fixed since the phenomenon studied here is quite localized to the impact area. Rock length (L r ) 5 mm

Fluid particle size 1 mm

In order to account for the centrifugal forces undergone by the Pelton bucket before the impact takes place, a static analysis of the bucket is carried out that precedes the impact transient analysis will involve the impact of the rock and the water jet on the Pelton bucket. Since the solid solver used for the impact calculations, Europlexus, is intended for transient dynamic calculations and uses the explicit Newmark central difference time scheme, no stiffness matrix or solver is available and hence no way to obtain the steady state response of the rotating bucket is possible. Thus, in order to obtain the displacement field of the bucket after reaching a steady state and just before the impact we carry out a static analysis using Code Aster, a finite-element solver developed by Electricité de France [EDF 00]. This solver, by default uses an implicit Newmark time-integration scheme for transient calculations and is capable of carrying out static analysis of struc-4. Numerical Simulations tures. In addition, the results can be written on a file format that can be subsequently used as an input file with Europlexus to launch the fluid-structure transient calculations that are of interest in the scope of the current thesis. Taking a diameter of D w = 0.40 m for the Pelton wheel and using the equations linking the tangential velocity of the bucket to the angular velocity of the turbine given in the Introduction, the angular velocity is set at ω = 500 rad/s. The displacement field of the static analysis are presented in Fig. 4.18. Once the steady state analysis is carried out, the fluid-structure transient analysis is launched using the displacement field of the static calculation as the initial configuration of the FSI calculation. The FSI calculation is launched using 2 computational cores for the fluid sub-domain, one core for the coupling software and one core for the solid solver.

Due to the complexity of the mesh of the solid sub-domain, the stability time-step for the solid sub-domain is fixed at ∆t s = 1.0 × 10 -8 s and the fluid one at ∆t f = 1.0 × 10 -7 s, hence a time-step ratio of m = 10 will be used for this calculation.

The results of the simulation are given in In Figure 4.19, the left column shows the case where the rock impacts the Pelton bucket in a dry configuration. The column on the right displays the fully coupled case. The indicator for damage in the steel bucket is once again the total equivalent plastic strain. In contrast to the case where a spheric projectile impacted a flat rectangular target head on, the damage caused by the projectile in the current coupled case is more pronounced than the damage caused by the dry impact. This phenomenon can be explained by carefully analysing the way the impact takes place in each case. For the dry case, due to the orientation of the rock at the beginning of the simulation, the actual surface of the the stone coming into contact with the Pelton runner is quite small. Besides, the contact time is rather short since, as shown in Figure 4.19, the stone immediately bounces off the surface of the bucket upon contact.

In contrast, for the fully coupled case, the hydrodynamic shape of the bucket prevents a thick water cushion to form between the two solid bodies and thus the speed of the stone is barely reduced before impact with the structure (v s ≈ -98 m/s). Additionally, the momentum of the water jet is transmitted to the stone and prevents it from quickly bouncing off from the bucket as was the case in the dry impact case. As such, the impact in the coupled case is extended both in time and space and thus the damage on the structure is greater than for the non-coupled case for this particular configuration.

Erosion damage can be visualized if elements surpassing 10% of cumulated plastic strain are considered as eroded. 4.20 shows the extent of the damage dealt after the impact in both the dry and the fully coupled case. As already explained, the extent of the damage is larger in the latter case due to a more persistent and larger contact area induced by the presence of the water jet which both guides the rock against the bucket and continues to transmit its momentum unto it.

For the FSI calculation presented here, the rock carried by the water jet was randomly placed near the extremity of the jet. Hence, the impact on the blade took place on a rather sturdy area of the bucket where the impact, although violent, would cause limited damage only. In order to consider a more unfavourable case scenario where the rock would impact a more fragile zone of the bucket, another simulation has been carried out where the stone is placed closer to the center of the jet. Due to the repositioning of the stone within the jet, the impact will initially take place on the converging blade that divides the jet into two halves known as the splitter. The results of this second simulation are presented in Fig. 4.21. Once again a comparison between the fully coupled case and a dry impact case are shown.

As one would expect, the damage induced on the target structure is much more pronounced this time around than in the previous case. In fact, by taking into account the of total equivalent plastic strain as indication of the damage on the turbine's bucket, an increase of around 50% for the maximum value is observed for the fully coupled case and there is around 100% increase for the maximum value in the dry impact case.

In contrast to the previous impact simulation, the dry impact case produced more damage on the structure than the fully coupled case. This effect can be explained if the way the impact takes place is analysed for each case. For the fully coupled scenario, there is a pressure spike that forms as the water jet encounters the splitter and is suddenly slowed down forming a water cushion that absorbs a bigger amount of the kinetic energy of the rock projectile. Additionally, in contrast to the previous simulation, due to the orientation of the splitter and the projectile upon impact, the surrounding water jet tends to drive the projectile away from the impact zone instead of thrusting it into the structure as was the case previously. As well, since the impact takes place before the water jet is split into two halves, the effect of the flux conservation is not present and hence the velocity of the water surrounding the jet is the same as the initial velocity of the jet. Hence, the impact velocity is v r ≈ 75 m/s before impacting the splitter this time around. In contrast, for the uncoupled simulation, the projectile hits the turbine on a weaker region at full speed coming closer to the case scenario presented in the previous sections where a projectile struck a flat plate head-on.

123

In order to better appreciate the extent of the damage that such an impact would induce on the bucket of the Pelton turbine, once again elements surpassing 10% damage will be considered eroded. Fig. 4.22 presents the erosion extent for both the dry impact case as well as the fully coupled case. As shown in Fig. 4.22, the extent of the erosion is larger in the case where no water jet is present in the calculations. However, for both cases, the extent of the damage is much more important with respect to the case where the impact took place in a more sturdy region of the bucket as was the case in the impact calculations presented earlier. By having a visual representation of the damage the impact of a stone can produce after such a violent impact without having to carry out costly experimental trails, the turbine manufacturer can determine if the damage induced on the turbine comes from the presence of large stones in the flow, which shouldn't be the case when operators follow the established security procedure.

Finally, the proposed multi-time step coupling technique allows to carry out fully coupled FSI simulations within a reasonable amount of time even for complex threedimensional cases as the one featured here. The present simulations were carried out with a time-step ratio of m = 10. In order to compare the duration of a simulation when using a multi time-step approach to a mono time-step simulation, the latter type of simulation has been carried out for the first 10% of the total duration of the simulation. Using this information, the computational time of a mono time-step simulation has been estimated.

The total computational times for each a mono time-step and a multi time-step calculations are presented in Table 4.11. The calculations presented here were launched using the same server as the one used for the impact calculations in 2-D. Two CPUs were used for the fluid solver, and one CPU each for the solid solver and the coupling program respectively. From Table 4.11, the use of the multi time-step coupling technique allows to obtain results for the coupled calculation around five times faster. Hence, the interest of using the proposed multi timestep coupling strategy is highlighted for complex 3-D calculations like the present one.

Conclusion

Validation examples for the heterogeneous time-step fluid structure coupling technique have been presented in the current chapter. The examples in 1-D show the basic properties of the proposed method such as its accuracy and convergence when smaller time-steps are used either in the fluid or in the solid. The 2-D examples focused on cases where the smaller time-step is used in the solid and highlighted how the proposed method can be applied even in cases of strong geometrical non-linearities. The validation examples allowed to compare the results of coupled FSI cases using the proposed method to the results obtained by other authors.

After the validation example phase, 2-D and 3-D examples of coupled calculations with impacts between a rock projectile modelled with quartz material properties and a steel target with the material properties of a Pelton bucket where carried out. For quantifying the damage induced on the steel target, the total equivalent plastic strain was observed. For 2-D and 3-D cases where the impact plane was perpendicular to the jet velocity vector carrying the projectile, the amount of damage induced on the target plate was much less pronounced when water jet was present as opposed to when the impact was done in a dry configuration. The reason for this was observed to be the influence of the water cushion formed between the projectile and the target, greatly absorbing much of the kinetic energy of the rock before impact. However, in the last simulation, where more realistic geometries where used, due to the configuration of the projectile and steel plate and to the hydrodynamic properties of the latter, the water cushion influence is less pronounced and the velocity of the projectile is almost the same as in the configuration where no jet is present. Additionally, the presence of the jet transfers its momentum and drags the projectile along the surface of the bucket, hence the damage observed on the latter is much more spread-out and pronounced in the coupled case.

General conclusion & perspectives

The work presented in the current thesis has focused in providing a tool to couple external specialized solvers in order to carry out transient fluid-structure interaction simulations.

The coupling of the explicit solid solver Europlexus, developped jointly by the French CEA and the European Commission's JRC, and the fluid solver Asphodel, developed by Andritz Hydro, was done by developing a coupling program handling the calculation of the pressure at the interface in order to correct the velocities at the boundary separating the fluid and the solid. The coupling of explicit solvers done in the current work adds a great deal of flexibility since the time integration scheme and the time-step size can be different and can be made compatible with calculations involving impacts between disjoint solid structures.

The first chapter of this thesis gave an overview on the different ways to solve fluidstructure interaction problems as well as on domain decomposition techniques and solid contact detection models which are at the base of the coupling technique used in the current work. From what was presented in this chapter, it was argued how the most convenient way to carry out the targeted FSI simulations was by coupling a mesh-less fluid solver, allowing to easily track the free surface and the fluid-solid interface to an explicit Finite-Element solid solver capable of handling contacts between solids. To overcome the time-step size constraints imposed by the explicit time integration schemes used by the solid and fluid solvers, inspiration from FETI-based techniques would be drawn in order to propose a multi-time step coupling strategy allowing to couple solvers using different time-steps.

The second chapter presented the governing equations for each the fluid and solid sub-domains and their discretization schemes in time and space. For the fluid, the basics behind the SPH-ALE technique were presented, which compared to the traditional SPH method, allows to impose a numerical particle transport velocity different than the actual flow velocity as well as better estimate intermediate states between particles are found by solving a one-dimensional Riemann problem, thus providing an appropriate framework to conveniently impose boundary conditions and normal velocity continuity in the case of fluid-solid coupled calculations. The widely used finite-element method was then briefly presented as the technique used for discretizing the solid governing equations in space.

Subsequently, the coupling between both discritized sub-domains was shown to be carried out by ensuring normal velocity continuity at the interface by using Lagrange multipliers. Regarding time integration, it was pointed out that the different time integration schemes are used for each sub-domain. Hence, for the fluid sub-domain, a explicit Runge-Kutta 2 (midpoint version) was applied, while the explicit Newmark Central Difference scheme was used in the solid sub-domain. Finally, the calculation of contact forces between solid meshes were deemed compatible with the FSI coupling tool and a brief description allowing to obtain them was also put forward.

The last two chapters of the current thesis focused on describing and validating the multi-time step fluid-structure interaction coupling strategy used throughout this work.

The third chapter introduced and described how the multi-time step coupling strategy can be used with smaller time steps in either the fluid or in the solid. The proposed coupling technique was based on existing solid sub-domain coupling techniques allowing to couple solid sub-domains using different time-steps and time integrators. Algorithms and flowcharts detailing the implementation of each of the coupling techniques were given and a stability proof was put forward in the case in which the smaller time-steps are found on the solid. On the last chapter, validation numerical examples were given comparing the results obtained with the proposed multi-step coupling technique to the results obtained by other authors or using the mono time-step coupling technique that was previously developed. Finally, coupled and uncoupled 2-D and 3-D problems aiming at quantifying the damage undergone by a steel target structure struck by a quartz rock projectile found within a water jet where carried out. The influence of the size of the jet as well as the geometrical form of the target and the projectile were determined to be important factors regarding the severity of the damage induced on the target structure after the impact took place.

The outlook for further research stemming from the current work is rather vast. This work shows that it is possible to couple fluid and solid sub-domains in highly non-linear cases using different time-steps and time integrators while being compatible with the simulation of contacts and impacts. Although some insight is given regarding the influence of the size of the projectile and the presence of a water jet on the damage undergone by the target structure, more simulations aiming at reproduce this phenomenon as close as possible to the real-life operating conditions of a Pelton turbine would be desirable. In order to do so, more ambitious simulations involving more than just a single bucket of the turbine could prove to be interesting. However, the increasing number of degrees of freedom in the solid and the fluid sub-domains would call for larger computational resources and development of algorithms capable of taking advantage of parallel computer architectures. The fluid solver is currently capable of running on parallel architectures, however, the solid solver, although capable of running in parallel as well, has not been rendered compatible with the coupling procedure in the case of its use for coupled fluid-structure calculations. Furthermore, the coupling software uses sparse matrices and direct solvers in to simplify and speed up the retrieval of the Lagrange multipliers. If larger algebraic systems are to be solved as a result of the larger numerical fluid and solid numerical models, the use of direct algorithms to solve the algebraic system becomes inefficient and the use of iterative, parallel capable solving algorithms becomes a necessity.

Regarding the validity of the models used to obtain a numerical solution of the coupled problems, as discussed over the course of this work, some simplifications have been considered in order to ease up the solution process. For example, the governing equations for the fluid-sub-domain, do not take into account viscosity-related terms and as a result, the fluid-solid coupling is only done using the normal components of the velocity. It would be thus interesting to develop a coupling algorithm allowing to take into account the tangential components of the velocity at the interface as this would allow to have some more insight into the effect of viscous induced forces on the outcome of the coupled calculations.

Considering the solid sub-domain, a material model capable of simulating sudden brittle fracture allowing the rock projectile to segment into tinier pieces would provide a way to better assess the extent of the damage induced by the impacting rock on the steel target structure. To carry out simulations with such a material model would be quite challenging since the surfaces that can become potentially wet after fracture would be changing continuously as the projectile sections into smaller pieces and thus the size of the interface problem can become quite large very quickly without knowing a priori the concerned degrees of freedom. A solid SPH model of the rock could also be used to ease the FSI treatment in the case of rock fragmentation. However a coupling strategy should be employed to link the domain discretized using the solid SPH to the domain discretized by the fluid SPH formulation.

Finally, erosion models allowing to predict damage on the structure after impact can be made compatible with the coupled calculations. For the proposed coupling technique, post-treatment techniques based on the total equivalent plastic strain can give an insight to the mesh zones that are prone to be eroded after an impact. For less violent phenomena in which no plastic state is reached or if it is negligible, an algorithm allowing to estimate damage through cyclic loads can be envisaged based on the results of the coupled calculation.

Introduction 1 11. 4 . 4 2. 3

 1443 Domain decomposition methods for fluid-structure interaction problems 9 1.1 Introduction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . 11 1.2 Coupling of fluid-structure solvers . . . . . . . . . . . . . . . . . . . . . 12 1.3 Fluid-structure interaction -the different approaches . . . . . . . . . . . . 15 1.3.1 Analytical Methods for fluid-structure interaction . . . . . . . . . 15 1.3.2 Semi-discretization methods for fluid-structure interaction . . . . 21 1.3.3 Coupling of semi-discretized methods for fluid-structure interaction 26 1.4 Domain decomposition methods -application to fluid-structure interaction problems .. . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . 31 1.4.1 Schur's primal method . . . . . . . . . . . . . . . . . . . . . . . 35 1.4.2 Schur's dual method . . . . . . . . . . . . . . . . . . . . . . . . 36 1.4.3 Multi-scale methods in time . . . . . . . . . . . . . . . . . . . . 37 Conclusion on the domain decomposition methods and its applicability to fluid-structure interaction problems . . . . . . . . . . . 39 1.5 Solid contact-impact modelling . . . . . . . . . . . . . . . . . . . . . . . 40 1.5.1 Impenetrability conditions . . . . . . . . . . . . . . . . . . . . . 41 1.5.2 Contact detection and contact force determination using the FEM 42 1.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 i List of Figures 1.14 The pinball algorithm uses circles or spheres created ad-hoc for contact detection while minimizing domain interpenetration . . . . . . . . . . . . 2.1 Allure of a kernel function used to weigh interactions between SPH particles 2.2 Local coordinate system to on which the solution to the Riemann problem will be sought . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Riemann problem solution space in the x ′ -t plane. . . . . . . . . . . . .

4. 1 4 . 4

 144 1-D bar coupled with a column of water -propagation of shock wave across the fluid-structure interface . . . . . . . . . . . . . . . . . . . . . 4.2 1-D propagation of shock wave across the interface -Comparison between the results obtained when integrating with same and different time-steps the two domains (smaller time-step in the fluid) . . . . . . . . . . . . . . 4.3 1-D propagation of shock wave across the interface -Comparison between the results obtained when integrating with same and different time-steps the two domains (smaller time-step in the solid) . . . . . . . . . . . . . . Kinematic continuity is prescribed at the micro-time scale when using the GC technique. Exact equilibrium is only prescribed at the macro timescale (black circles), an interpolation allows to prescribe equilibrium at the micro time-scale (red crosses) . . . . . . . . . . . . . . . . . . . . . 4.5 2-D hydrostatic water interacting with a linear elastic plate . . . . . . . .

Figure 1 :

 1 Figure 1: Erosion damage visible on a Pelton turbine bucket at the end of its functional lifetime

Contents 1. 1 1. 3 . 3

 133 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . 1.2 Coupling of fluid-structure solvers . . . . . . . . . . . . . . . . . . . . . 1.3 Fluid-structure interaction -the different approaches . . . . . . . . . . 1.3.1 Analytical Methods for fluid-structure interaction . . . . . . . . . . 1.3.2 Semi-discretization methods for fluid-structure interaction . . . . . Coupling of semi-discretized methods for fluid-structure interaction 1.4 Domain decomposition methods -application to fluid-structure interaction problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4.1 Schur's primal method . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1 . 3 :

 13 Figure 1.3: Cylinder within an unbounded potential flow

a) M = ρ 0 FE

 0 /ρ 0 S known as the mass number since it is the ratio of densities of the fluid and the solid b) C Y = ρ 0 F U 2 0 known as the Cauchy number which measures the magnitude of dynamic pressure-induced deformations.

Figure 1 . 4 :

 14 Figure 1.4: Illustration of two case scenarios in which dimensional analysis can be used to solve the fluid-structure interaction problems. A solid floats over a fluid at rest (a). An airfoil inside an airflow (b)

FFigure 1 . 6 :

 16 Figure 1.6: Domain Ω s discritized with a FEM mesh under different types of boundary conditions

Figure 1 . 7 :

 17 Figure 1.7: FVM cell Ω f v and the associated fluxes along the cell walls are in

Figure 1 . 8 :

 18 Figure 1.8: SPH particle i, fluxes are integrated over the kernel support of the particle

Finite

  volume and finite-element Eulerian meshes can be used as well to model free surfaces in fluid flows and to keep track of the fluid-solid interface for fluid-structure interaction simulations. In order to track the fluid interfaces, techniques such as the Volume of Fluid (VOF) [HIR 81] or Level-Set [OSH 06] can be used. The former technique keeps track of the volume fraction of fluid inside the different cells of the mesh in order to locate the position of the free surface in question. The Level-Set technique replaces the volume fractions by using a distance to the interface function. In [HAC 15], a monolithic fluidstructure interaction technique is presented where the level-set method is used along with a dynamic adaptive mesh algorithm in order to enhance the detection of discontinuities at the fluid-solid boundary.

U 0 s

 0 Figure 1.9: NACA airfoil within air flow. ALE fluid mesh can follow solid motion

Figure 1 . 10 :

 110 Figure 1.10: Decomposition possibilities for a domain Ω into smaller sub-domains Ω 1and Ω 2 separated by interface Γ 12 . On the right decomposition, the interface's dimension is the same as the dimension of Ω. On the decomposition in the left, the dimension of the interface is decreased by one

F

  d are the boundary forces applied on ∂ F Ω, f d are the body forces applied on Ω and U d are the kinematic boundary conditions applied on ∂ u Ω.

  1. Domain decomposition methods for fluid-structure interaction problems pliers. It is clear how the interest behind the domain decomposition techniques presented in the previous section is amplified if these can be extended to the case of transient analysis where each sub-domain involved can be potentially integrated with a different time-scheme and time-step. Research in this sense has sprouted a vast domain where the interface continuity equations can either be satisfied at the smaller of time-scales, as put forward mainly by the GC method [GRA 01], or at the larger time-scale, as seen in the methods proposed by the PH method[PRA 04] and the MGC method [MAH 11]. More recently, Brun et al. [BRU 15] proposed a general framework for building energy conserving methods at the two time scales (kinematic continuity at the large or the fine timesteps), by cancelling the interface pseudo-energy involved in the so-called energy method, this was originally proposed by Hughes for proving the stability of hybrid explicit-implicit time integrators [HUG 12].

Figure 1

 1 Figure 1.12: Solid domains Ω 1 and Ω 2 in contact over surface Γ c

Figure 1 . 13 :

 113 Figure 1.13: The sliding lines algorithm allows to determine contact using the mesh elements and nodes

  Figure 1.14: The pinball algorithm uses circles or spheres created ad-hoc for contact detection while minimizing domain interpenetration

Figure 2 . 1 :

 21 Figure 2.1: Allure of a kernel function used to weigh interactions between SPH particles

Figure 2 . 2 :

 22 Figure 2.2: Local coordinate system to on which the solution to the Riemann problem will be sought

Figure 2 . 3

 23 Figure 2.3 shows the solution space for the Riemann problem is divided into four main regions: the left region W L , the right region W R , the left star region W * L and the right star region W * R . The wave associated with the λ 2 value can be considered a contact discontinuity. along which the densities and the velocities in the normal directions have the same values (ρ * L= ρ * R = ρ * and v * L,x ′ = v * R,x ′ = v * x ′ ).The remaining two waves associated to eigenvalues λ 1 and λ 3 can either be a shock-wave or a rarefaction wave. The values of the conservation variables are known for W L and W R since they

Figure 2 . 3 :

 23 Figure 2.3: Riemann problem solution space in the x ′ -t plane.

  problem. In 2007, Marongiu proposed an original way to tackle the way boundary conditions are imposed by writing them as a solution to a partial Riemann problem between a solid boundary and the fluid particles whose kernel is sectioned by the former [MAR 07].

Figure 2 . 4 :

 24 Figure 2.4: Reference problem for establishing the solid sub-domain's governing equations

Figure 2 . 5 :

 25 Figure 2.5: Illustration of the outward pointing normals, n s and n f , for each sub-domain

  52) where α = 1, 2 denotes a principal direction in the 2-D Cartesian plane and I is one of the DOF of the problem, I ∈ (1, N N ). Hence, f ext f (I) is the value of the I th DOF in the f ext f 2. Fluid-solid coupling strategy vector.

Contents 3 . 1

 31 Incompatible time-step implementation . . . . . . . . . . . . . . . . . . 78 3.1.1 The GC method with smaller time-steps in the solid and large interface displacements . . . . . . . . . . . . . . . . . . . . . . . . . 79 3.1.2 Smaller time-steps in the fluid . . . . . . . . . . . . . . . . . . . . 87 3.2 Stability of the algorithm. . . . . . . . . . . . . . . . . . . . . . . . . . 93 3.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . 96

3 .

 3 FSI Coupling with different time scales and solid impacts tial integration [MAR 07]: the resulting dissipation is usually large. Lastly, the operators needed for the interface calculation are much easier to obtain when compared to other multi-time-step domain methods[MAH 11].

Figure 3 . 1 :

 31 Figure 3.1: Overview of the coupling procedure when the solid uses a smaller time-step than the fluid, i.e. ∆t f = m∆t s

Figure 3 . 2 :

 32 Figure 3.2: Flowchart broadly describing how the algorithm is implemented for the case when smaller time-steps are used in the solid sub-domain

  H j and b j to obtain Λ j and sends it to SLS -SLS calculates a

Figure 3 . 3 :

 33 Figure 3.3: Overview of the coupling procedure when the fluid uses a smaller timestep than the fluid, i.e. ∆t s = m∆t f and information from the interface comes only at the beginning of each macro-time-step t 0 = t n

Figure 3 . 4 :

 34 Figure 3.4: Flowchart broadly describing how the algorithm is implemented for the case when using smaller time-steps in the fluid sub-domain

  CPL and SLS interact to obtain U

  42) Using equation (3.18) we can rewrite equation (3.42) as:

4. 1

 1 .1 1-D propagation of shock wave across the fluid-structure interface For the first test case we couple cantilever 1-D linear bar to a water tube inside which a strong pressure gradient induces a shock wave across the interface (Figure 4.1).

Figure 4 Figure 4

 44 Figure 4.1: 1-D bar coupled with a column of water -propagation of shock wave across the fluid-structure interface

Figure ( 4 .

 4 Figure(4.2) shows the comparison between the results of the simulation using the method described in the current paper and another one in which both sub-domains are

Figure 4 . 3

 43 Figure 4.3 shows the result comparison between the proposed method and those obtained from a uniform time-step simulation carried out by [ZHE 13] (∆t s = ∆t f = 10 -8 ).Once again, we see that good agreement also exists when using a smaller time-step to integrate the solid through the use of the GC method for time-steps that differ considerably.

Figure 4 . 3 :

 43 Figure 4.3: 1-D propagation of shock wave across the interface -Comparison between the results obtained when integrating with same and different time-steps the two domains (smaller time-step in the solid)

Figure 4 . 4 :

 44 Figure 4.4: Kinematic continuity is prescribed at the micro-time scale when using the GC technique. Exact equilibrium is only prescribed at the macro time-scale (black circles), an interpolation allows to prescribe equilibrium at the micro time-scale (red crosses)

  g

Figure 4 .

 4 Figure 4.5: 2-D hydrostatic water interacting with a linear elastic plate

  top-central node Δt f = 10Δt s same time-step

Figure 4 .

 4 Figure 4.6: 2-D hydrostatic water on a linear elastic plate -Comparison between the results obtained when integrating with same and different time-steps in both domains

Figure 4 . 7 :

 47 Figure 4.7: Initial configuration of the test case: breaking dam flow on an elastic wall

  .12 0.14 0.16 0.18 0.2 0.22 0.24 0.26 0.28 0.3 Displacement (m) time (s) Horizontal displacement of the top-left node Europlexus mono time-step Europlexus multi-time step Code Aster mono time-step Walhorn et al. Idelsohn et al.

Figure 4 . 9 :

 49 Figure 4.9: Result comparison of the displacement of the top left node of the structure until t = 0.3 s

Table 4 . 5 :- 4 m

 454 Geometric and discretization 2-D jet and plate Water width (L f ) 0Plate mesh spacing 5 × 10 -4 m

Table 4 . 6 :

 46 Simulation parameters for the mesh convergence analysis Domain ∆x 1 ∆x 2 ∆x 3 ∆t fluid 0.5 mm 0.25 mm 0.1 mm 5 × 10 -8 s solid 0.5 mm 0.25 mm 0.1 mm 1 × 10 -8 s

Figure 4 .Figure 4 .

 44 Figure 4.11: 2-D jet impact on plate, mesh convergence study

  ine the extent of the damage on the structure caused by ductile fracture [KIM 15, SEI 16]. The material parameters for the steel plates are taken from Andritz official documentation.

Figure 4 .

 4 Figure 4.12 presents the geometrical configuration of the simulation and Table4.7 presents the geometrical parameters taken into consideration for this simulation.

Figure 4 . 4 m

 44 Figure 4.12: 2-D jet and stone impact on plate

Finally, the

  fluid parameters include its reference density ρ 0 f = 1000 kg/m 3 and its reference speed of sound c 0 f = 1000 m/s. The results of the simulation are presented in the following Figures.

Figure 4

 4 Figure 4.13: Impact of a quartz rock on a steel plate. Comparison of the damage induced through a dry impact and through a fully coupled impact

. 13 :

 13 Figure 4.13: Impact of a quartz rock on a steel plate. Comparison of the damage induced through a dry impact and through a fully coupled impact

Figure 4

 4 Figure 4.13 presents a comparative result of the damage undergone by the steel plate both in the case when the rock impact occurs in a "dry" configuration (top row) as well as in a fully coupled configuration (bottom row). The results presented in this Figure

Figure 4 .

 4 Figure 4.14: Impact of a quartz rock on a steel plate. Influence on the size of the stone on the damage undergone by the steel plate.

Figure 4 .

 4 Figure 4.15: 2-D jet and stone impact on plate, comparison between mono time-step and multi time-step calculations

Figure 4 . 16 :

 416 Figure 4.16: Impact of a quartz rock on a steel plate. Comparison of the damage induced through a dry impact and through a fully coupled impact

Figure 4 .

 4 Figure 4.17: Pelton turbine bucket top and side views

Figure 4 . 18 :

 418 Figure 4.18: Steady state displacement field of the Pelton Bucket due to the centrifugal forces induced by its rotation

Figure 4

 4

Figure 4 .

 4 Figure 4.19: Impact of rock on Pelton turbine bucket

Figure 4 .

 4 20 presents the comparison of the surfaces eroded for both the dry and the fully-coupled case.

Figure 4 .

 4 Figure 4.20: Erosion damage on the Pelton turbine bucket

Figure

  Figure4.20 shows the extent of the damage dealt after the impact in both the dry and the fully coupled case. As already explained, the extent of the damage is larger in the latter case due to a more persistent and larger contact area induced by the presence of

Figure 4 .

 4 Figure 4.22: Erosion damage on the splitter section of the Pelton bucket

Table 4 . 11 :

 411 Elapsed calculation time for same time-step and multi time-step coupling calculations for the 3-D impact problem time-step ratio calculation time m = 1 43 hrs (estimated) m = 10 9 hrs

  FOLIO ADMINISTRATIF THESE DE L'UNIVERSITE DE LYON OPEREE AU SEIN DE L'INSA LYON NOM : NUNEZ RAMIREZ DATE de SOUTENANCE : 29/05/2017 Prénoms : Jorge Israel TITRE : A multi time-step partitioned approach for the coupling of SPH and FE methods for nonlinear FSI problem NATURE : Doctorat Numéro d'ordre : 2017LYSEI040 Ecole doctorale : Mécanique, Energétique, Génie Civil, Acoustique Spécialité : Mécanique RESUME : Dans le cadre de ce travail, une technique non-intrusive est proposée pour coupler la méthode Smoothed Particle Hydrodynamics (SPH) à la méthode des Eléments Finis (MEF) afin de résoudre numériquement des problèmes dynamiques et non-linéaires d'interaction fluide-structure en permettant l'utilisation des pas de temps différents dans les deux domaines de calcul (fluide et solide). Ces développements furent motivées par le besoin de simuler numériquement des phénomènes rapides et très non-linéaires qui prennent en compte des impacts en se servant des intégrateurs temporels explicites dans chaque sous-domaine de calcul (Newmark explicite pour le solide et Runge-Kutta 2 pour le fluide). De ce fait, le pas de temps de stabilité est limité par des caractéristiques intrinsèques au modèle numérique du phénomène étudié et en conséquence, il devient important de pouvoir intégrer chaque sous-domaine numérique avec un pas de temps égale ou proche de son pas de temps de stabilité. Pour permettre d'utiliser un pas de temps proche du pas de temps de stabilité pour chaque sous-domaine, des méthodes de décomposition de domaines dual-Schur sont implémentées et validées pour des cas en 1-D, 2-D, et 3-D. Des simulations numériques des impacts des cailloux sur des aubes des turbines hydrauliques sont aussi effectuées afin de prédire le dommage que cet évènement peut occasionner sur cette dernière pendant son fonctionnement. La méthode de couplage fluide-structure proposée par ce travail est utilisée à la fois par trois autres travaux doctoraux pour prédire le dommage sur des aubes des turbines hydrauliques résultant de l'impact répété des gouttes d'eau, l'érosion par cavitation ainsi que l'érosion hydro-abrasive provoquée par les petits sédiments présents dans l'´ecoulement. MOTS-CLÉS : COUPLAGE FLUIDE-STRUCTURE, ELEMENTS FINIS, SMOOTHED PARTICLE HYDRO-DYNAMICS, IMPACTS, DECOMPOSITION DES DOMAINES Laboratoire (s) de recherche : LaMCoS Laboratoire de Mécanique des Contacts et des Structures Directeur de thèse: Alain COMBESCURE Président de jury : Composition du jury : Mhamed SOULI, Elie HACHEM, Olivier ALLIX, Thouraya BARANGER, Jean-christophe MARONGIU, Michael BRUN, Alain COMBESCURE

  and CPL send real I s status to FLS

	else
	-SLS directly obtains U m s , sends I s to FLS
	end if
	end if
	end for
	end for

Table 4 . 1 :

 41 Computed order of accuracy in time for the mono time-step algorithm and the multi time-step algorithm with smaller time-steps in the solid

	time-step ratio	τ	τ 1	τ 2
	m = 1	Ω(u s (L 0 s )) 1.949 2.044
	m = 100	Ω(u s (L 0 s )) 1.605 1.871

Table 4 . 2 :

 42 Compared observed order of accuracy in time for the mono time-step algorithm and the multi time-step algorithm with smaller time-steps in the fluid

	time-step ratio	τ	τ 2

Table 4 . 3 :

 43 Simulation parameter used in the 2-D hydrostatic test case.

	Water width (L)	1 m
	Water height (H)	2 m
	Elastic wall thickness (e)	0.05 m
	Fluid particle size	0.01 m
	Rigid wall spacing	0.01 m
	Elastic wall mesh spacing	0.01 m
	Fluid reference density	1000 kg/m 3
	Fluid reference sound speed	50 m/s
	Solid Young's modulus	67.5 GPa
	Solid density	2700 kg/m 3
	Poisson coefficient	0.34
	Solid 2-D hypothesis	plane strain

Table 4 .

 4

	8: Elapsed calculation time for same time-step and multi time-step coupling
	calculations
	time-step ratio calculation time
	m = 1	31 min
	m = 10	12 min

Table 4 .

 4 9: Geometric and discretization parameters for the 3-D impact problem

	Water width (L f )	0.1 m
	Plate width (L s )	0.2 m
	Plate height (h s )	0.05 m
	Rock radius (r s )	20 mm
	Fluid particle size	5 × 10 -3 m

Cette thèse est accessible à l'adresse : http://theses.insa-lyon.fr/publication/2017LYSEI040/these.pdf © [J. Nunez Ramirez], [2017], INSA Lyon, tous droits réservés

Multiphysics coupling techniques aimed at fluid-structure interaction applications

Domain decomposition techniques aimed at the coupling of sub-domains with heterogeneous time integrators and different time scalesIn general, for all of the coupling and domain decomposition methods that will be presented in this section, when splitting a numerical domain Ω into, for example, two smaller domains Ω 1 and Ω 2 sharing a common interface Γ, in order to obtain a perfect link between the domains, one must ensure for the kinematic variable u i Γ and the dynamic variable F i Γ the following:(i) Kinematic continuity between sub-domains :u 1 Γ = u 2 Γ

Cette thèse est accessible à l'adresse : http://theses.insa-lyon.fr/publication/2017LYSEI040/these.pdf © [J. Nunez Ramirez], [2017], INSA Lyon, tous droits réservés Domain decomposition methods -application to fluid-structure interaction problems

Cette thèse est accessible à l'adresse : http://theses.insa-lyon.fr/publication/2017LYSEI040/these.pdf © [J. Nunez Ramirez], [2017], INSA Lyon, tous droits réservés 1. Domain decomposition methods for fluid-structure interaction problems

Cette thèse est accessible à l'adresse : http://theses.insa-lyon.fr/publication/2017LYSEI040/these.pdf © [J. Nunez Ramirez], [2017], INSA Lyon, tous droits réservés 2.5 Conclusion on FSI couping with the SPH-FEM methods . . . . . . . . . 74

Remerciements

Chapter 4 Numerical Simulations

In this chapter, validation examples of the multi-step coupling strategy will be given. The examples will feature 1-D, 2-D and 3-D cases, where the accuracy of the proposed method will be assessed and, for some cases, comparisons with the results of other authors will be established. Following the validation examples, full simulations of the impacts of stones on turbine blades will be featured showing the coupling strategy's capacity to carry out such complex simulations and giving an insight to the extent of the damage such as phenomenon can induce on a turbine blade.