Renée El

Melhem Bat

Blaise Pascal

Insa : R Gourdon

M Stéphane

M Gérard Scorletti
email: gerard.scorletti@ec-lyon.fr

Keywords: Érosion, Couplage fluide-structure, Fatigue, Éléments Finis, Smoothed Particle Hydrodynamics Droplet impact, Erosion, Fluid-structure interaction, Fatigue, Finite Element Method, Smoothed Particle Hydrodynamics

The goal of this work is the comprehension and the numerical simulation of erosion caused by repeated droplet impact on Pelton turbine buckets, to predict the lifetime of these components. First, waves propagation phenomenon inside fluid and solid domains are presented, which allows determining the time evolution and spatial distribution of contact pressure, and the birth of lateral microjets by supersonic ejection of the fluid on the contact. Experimental studies of erosion by droplet impact highlight a fatigue crackingbased erosion mechanism. Then, coupled FSI computation are performed. The solid subdomain is discretized by the Finite Element Method (FEM), and the fluid subdomain by the Smoothed Particle Hydrodynamics (SPH), which is a particle method (meshless) effectively recommended for large distortions and free surface tracking. Stress analysis confirms the cyclic nature of the damage mechanism, and erosion simulation is performed using multiaxial fatigue criteria. The first selected criterion is a general one from the American Society of Mechanical Engineers (ASME) using principal values of stress differences over time. The second one is the Dang van 2 criterion, belonging to the family of critical plane criteria. This criterion considers separately the effects due to hydrostatic stress on one hand, and the ones induced by maximum local shear on the other. These two criteria are used to building the equivalent eroded zones of the solid subdomain for a given number of impacts, which allows to qualify this procedure as a predictive predictive. Finally, a parametric study for different droplet sizes and velocity is computed, and the effects of a coating layer are investigated.

Je tiens à adresser mes plus chaleureux remerciements à M. Alain Combescure, mon directeur de thèse, qui s'est montré attentif, bienveillant, compréhensif, tolérant, attentionné, patient, encourageant, extrêmement réactif et toujours disponible. Grâce à ses précieuses connaissances, son expérience et son recul, il m'a judicieusement aiguillé dans ma progression tout au long de cette aventure longue de trois années. En effet, chaque réunion de suivi s'est révélée stimulante et me motivait dans la poursuite du travail, même en étant bloqué, puisque les discussions menaient à des solutions claires et des démarches efficaces, le tout me permettant de m'approprier le sujet. Je considère que l'avoir eu comme Maître fut un immense privilège et une chance inouïe.

Je remercie particulièrement M. Jean-Christophe Marongiu, chef de l'équipe Turbine Physics chez ANDRITZ Hydro, pour avoir proposé et mis en place ce projet et nous avoir encadrés, pour sa disponibilité ainsi que ses explications claires concernant la dynamique des fluides numérique (notamment la méthode SPH) et le temps passé à nous aider à maitriser les outils numériques nécessaires.

Je remercie en outre les autres doctorants du projet PrEDHyMa, avec qui j'ai eu de nombreux échanges fructueux et scientifiquement pertinents : Wiebke Boden, Saira Pineda et tout particulièrement Jorge Nuñez, véritable ami, d'une grande aide, qui a généreusement partagé ses connaissances dans l'utilisation du shell sous Linux et des langages C et ForTran.

Je souhaite aussi exprimer ma gratitude envers M. Etienne Parkinson, Technical Manager et chef du département R&D chez ANDRITZ Hydro Vevey qui m'a permis de travailler dans un environnement favorable à l'accomplissement de ce travail de recherche, et aussi envers les autres collaborateurs du département R&D d'ANDRITZ Vevey et Villeurbanne qui m'ont chaleureusement accueilli lors de cette aventure : les deux Nicolas, Steven, Adrien, Hélène, Matthias, Martin, Magdalena, Ralph, ansi que Hervé pour son aide et ses judicieux conseils en programmation (particulièrement en C) ainsi que sur CMake.

Mes remerciements vont également aux doctorants et post-doctorants de l'ancienne équipe MSE du LaMCoS pour leur accueil lors de mes arrivées (en tant que docrorant "voyageur") : Jérémy, Espoir, Eric, Pipo, Damien, Emmanuel, Alexandre, Nacer, Nicolas, Yannick, Fatima, Guillaume, Salvatore, Corentin, Wenfeng, Yancheng, Kevin, Marie, Lv, ainsi qu'à ceux des nouvelles équipes MIMESIS et MULTIMAP : Tristan, Thibaut, Pierre, Matthieu, Arthur, Thomas, Zhi Kang, et tous ceux des autres salles et étages... Enfin, une grande part de ma réussite revient à mes parents qui m'ont encouragé et soutenu depuis toujours dans mes réalisations, et particulièrement durant cette "odyssée". 6

Résumé

Le but du travail présenté est la compréhension puis la simulation numérique des mécanismes d'érosion des augets de turbine Pelton par impacts répétés de gouttes d'eau dans le but de prédire la durée de vie des composants. Tout d'abord, les phénomènes de propagation d'ondes dans les milieux fluide et solide sont étudiés. Cela permet de mettre en lumière l'évolution temporelle et la distribution spatiale des pressions de contact, et l'apparition de microjets par éjection supersonique du fluide au contact. Les études expérimentales de l'érosion par gouttes d'eau traduisent un dommage basé sur la fissuration par fatigue. Des simulations numériques en dynamique rapide couplées fluide-structure sont alors effectuées. Le domaine solide est discrétisé par la Méthode des éléments Finis (MEF), et le domaine fluide par la méthode Smoothed Particle Hydrodynamics (SPH), qui est une méthode particulaire (sans maillage) particulièrement adaptée aux grandes distorsions et au suivi des surfaces libres. L'analyse des états de contraintes vient corroborer la nature cyclique de l'endommagement. La simulation d'érosion est alors réalisée à l'aide de critères de fatigue multiaxiaux. Le choix se porte vers un premier critère général de l'American Society of Mechanical Engineers (ASME), utilisant les valeurs principales des différences de contraintes au cours du temps. Le second choix concerne un critère à plan critique : le critère de Dang Van 2. Il traite séparément la contrainte hydrostatique et le cisaillement alterné maximal local. Ces critères permettent de définir les régions érodées du solide au bout d'un nombre d'impact donné, ce qui fait de cette démarche une méthode prédictive. Une étude paramétrique pour différentes tailles de gouttes et vitesses d'impact est ensuite réalisée, puis on évalue l'influence de la présence d'une couche de coating. Dang Van 2 criterion on right, and pressure inside fluid at t = 460 ns after impact for the case (φ = 1 mm,V = 100 m.s -1 ) . . . . . . . . . . . . . . 115 4.4 Number of cycles to failure according to ASME criterion on left and Dang Van 2 criterion on right, at t = 460 ns after impact for the case (φ = 1 mm,V = 100 m.s -1 ) . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 4.5 Influence of droplet diameter and impact velocity on maximum stress amplitude according to both criteria . . . . . . . . . . . . . . . . . . . . . . 117 4.6 Influence of droplet diameter and impact velocity on minimum number of cycles to failure according to both criteria . . . . . . . . . . . . . . . . . 117 4.7 Influence of droplet diameter and impact velocity on number of eroded elements according to both criteria for N lim = 5 • 10 7 cycles . . . . . . . . 118 4.8 Influence of droplet diameter and impact velocity on erosion depth according to both criteria for N lim = 5 • 10 7 cycles . . . . . . . . . . . . . . . . 119 4.9 Fatigue function and non-eroded elements for N lim = 5 • 10 7 cycles according to both criteria, at t = 460 ns after impact for the case (φ = 0.5 mm,V = 100 m.s -1 ) . . . . . . . . . . . . . . . . . . . . . . . . . . 120 vii 4.10 Fatigue function and non-eroded elements for N lim = 5 • 10 7 cycles according to both criteria, at t = 460 ns after impact for the case (φ = 0.5 mm,V = 200 m.s -1 ) . . . . . . . . . . . . . . . . . . . . . . . . . . 4.11 Fatigue function and non-eroded elements for N lim = 5 • 10 7 cycles according to both criteria, at t = 460 ns after impact for the case (φ = 1 mm,V = 100 m.s -1 ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.12 Fatigue function and non-eroded elements for N lim = 5 • 10 7 cycles according to both criteria, at t = 460 ns after impact for the case (φ = 1 mm,V = 200 m.s -1 ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.13 Fatigue function and non-eroded elements for N lim = 5 • 10 7 cycles according to both criteria, at t = 460 ns after impact for the case (φ = 2 mm,V = 100 m.s -1 ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.14 Fatigue function and non-eroded elements for N lim = 5 • 10 7 cycles according to both criteria, at t = 460 ns after impact for the case (φ = 2 mm,V = 200 m.s -1 ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.15 Dang Van 2 stress amplitude for standard mesh and refined one in the highest stress region, at t = 460 ns after impact for (φ = 1 mm,V = 100 m.s -1 )125 4.16 Dang Van 2 fatigue function and non-eroded elements for N lim = 10 7 cycles using standard mesh and refined one in eroded region, at t = 460 ns after impact for (φ = 1 mm,V = 100 m.s -1 ) . . . . . . . . . . . . . . . . 4.17 Comparison of hydrostatic stress between the pure basic material and the coated one at different frames for the case (φ = 1 mm,V = 100 m.s -1 ) . . 4.18 Comparison of Von Mises stress between the pure basic material and the coated one at different frames for the case (φ = 1 mm,V = 100 m.s -1 ) . . 4.19 Comparison of maximum value of stress amplitude over the mesh according to both criteria between the two material configurations for the case (φ = 1 mm,V = 100 m.s -1 ) . . . . . . . . . . . . . . . . . . . . . . . . 4.20 Comparison of minimum value of number of cycles to failure over the mesh according to both criteria between the two material configurations for the case (φ = 1 mm,V = 100 m.s -1 ) . . . . . . . . . . . . . . . . . . 4.21 Comparison of ASME stress amplitude between pure basic material and coated one at t = 460 ns after impact for φ = 1 mm and V = 100 m.s -1 . 4.22 Comparison of Dang Van 2 stress amplitude between pure basic material and coated one at t = 460 ns after impact for φ = 1 mm and V = 100 m.s -1 4.23 Comparison of number of cycles to failure according to ASME criterion between pure basic material and coated one at t = 460 ns after impact for φ = 1 mm and V = 100 m.s -1 . . . . . . . . . . . . . . . . . . . . . . . 4.24 Comparison of number of cycles to failure according to Dang Van

List of Figures

Nomenclature

The complete list of notations is given in Appendix A.

Subscripts and superscripts

Symbol Description

• i i line component of vector

• i j i line and j column component of matrix

Introduction

Hydraulic turbines can undergo severe damage during operation, because of low quality water or detrimental flow conditions. Damage induces maintenance costs and power production losses, and can also endanger safety of installations. Hydropower plants operators and turbine manufacturers are interested in extending overhaul periods by reducing damage intensity and protecting turbine components with surface treatments, but accurate and reliable prediction of damage is however missing. The present work is related to the erosion arising from repeated impacts of high speed water droplets on specific parts of Pelton turbines. Indeed for high head Pelton units, the jet of water is composed of a liquid core surrounded by droplets. Observations show that regions of impact of these droplets exhibit specific erosion patterns. The aim of the present work is to understand the mechanism of erosion by water droplet impact and build a numerical tool in order to simulate the corresponding damage, which allows predicting the lifetime of Pelton buckets.

This work is funded by the framework of the European project PrEDHyMa, a Marie Skłodowska Curie action fostering the collaboration between two French research institutions : Institut National des Sciences Appliquées de Lyon and École Centrale de Lyon, and a private Austrian Swiss-based corporation : ANDRITZ Hydro.

Hydroelectric power generation

Today, hydroelectricity appears to be one of the main sources for energy production of renewable nature. Currently, among all renewable energy sources, energy produced with hydroelctric stations represents 10% of energy generation and almost 17% of the global energy produced [ADI 16]. The exploitable sources for production of hydroelectricity have reach their peak for industrialized countries, this market is currently rising in developing countries. For the year 2015, the total energy output coming from hydroelctricity was estimated at 3940 TWh [ADI 16].

The main benefits of hydroelectricity making it an attractive and viable choice are the following :

• Its flexibility which comes from the ease and speed with which the power station can be stopped and restarted based on power need.

• The relative low cost of the energy produced because hydroelectric plants can be Introduction relatively inexpensive to build, they require few operative staff, and their expected functional lifetime can last for decades.

• The ability to stock the energy producing resources for use during peak demand periods. Moreover, these water reservoirs built in addition to the stations can dynamize the economic activity of the region such aquaculture or tourism.

Despite the advantages of hydroelectricity listed below, one must also consider some drawbacks :

• The construction of the power station leads to disruption to the local natural environment. Indeed hydroelectric station requires the flooding of the region to build a dam, which disturbs the local ecosystem and habitants must be displaced.

• The accumulation of silt in the premises of station and reservoir, transported by the flow of rivers and streams can saturate the reservoir which causes the the station to become non-operational.

The Pelton wheel

Hydroelectric power stations can produce energy with three main types of turbines : Francis, Kaplan, Pelton. The presented study concerns the last class : the Pelton turbine. They have been used since their inception by the American inventor and entrepreneur Lester Allan Pelton in the late 19 th century [WIL 74]. The main difference between Pelton turbine and centuries old water wheel is the number of blades/buckets in contact with moving water. Indeed, where the kinetic energy was transferred to the partially submerged blades only in old systems, the Pelton turbine allows to transfer the energy on each of the buckets at high velocity through a converging nozzle (Figure 1). The design of the bucket has also evolved in order to maximize the energy transfer and minimize losses. The contact of the water on bucket forms a U-turn and exits the bucket at much lower velocity, thus Pelton-based generators are single stage. Pelton turbines are commonly preferred for high dynamic heads, greater than 100 m, and low flow rates, where discharges are lower than 50 m 3 .s -1 . They come in a large variety of sizes : from large ones capable of outputting 400 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 V t is usually set to be the half of the jet velocity impacting the bucket [WIL 74] :

V t = È 2gH p 2
where g denotes the acceleration due to gravity field, and H p is the dynamic head. For a wheel with a diameter φ p , the angular velocity Ω r of the turbine runner is : For the current study, a high dynamic head will be considered as H p = 2000 m. Hence, the corresponding water jet velocity is V p = È 2gH p = 200 m.s -1 , and the tangential velocity of the runner will be V t = 100 m.s -1 . Thus, the relative velocity of the jet on the runner is V j→p = 100 m.s -1 , assuming the jet and tangential vector on contact point being collinear. The water droplets will have a velocity of the same order of magnitude as the jet and their diameter will be in the millimeter range, based on observations.

The PrEDHyMa project

The PrEDHyMa project, whose acronym stands for Prediction of Erosion Damage on Hydraulic Machines, consists in studying four possible sources of damage that hydroelectric turbines undergo during their lifetimes. The project focuses on buckets of Pelton runners, because these turbines are particularly used where high heads are available, notably in mountainous regions (Figure 2). The four sources studied in this project are :

• Gravels and stones impacts : The impact of gravels or stones found in streams during the turbine operation can lead to high damage on the buckets. This event should usually not occur, because precautions are normally taken to prevent it, but simulations of this kind of damage are important since the intensity of the impact can heavily damage the turbine.

• Droplets impacts : The impingement of water droplets on the surface of blades/buckets of an hydroelectric turbine can generate damage after repeated exposure time of this phenomenon. This kind of damage is different of the previous one, because erosion of the surface does not appear suddenly, but after prolonged exposure time.

• Hydro-abrasive erosion : The presence of small sediments in streams and rivers provides a source of erosion which is often responsible for the damaging of Pelton turbines. The damage is visible after prolonged exposure because of the low size of the solid sediments.

• Cavitation erosion : This type of damage is commonplace among various types of turbines. Cavitation phenomenon is produced by pressure drop inside the liquid core interacting with blades. Gas pockets form, then implode, creating a high-speed microjet that can damage turbines blades in the long run.

The current work deals with the second source of damage : droplet impact erosion.

Droplet impact erosion

An old research topic Moreover, the previous investigations about this erosion mechanism are built from experimental observations and the generated models are phenomenological ones. Erosion mechanism is described, but the physical phenomena responsible for it are not quantified and cannot even be estimated if the experimental tests are not available. Other studies have been performed with numerical simulations of the water droplet impacting a target. This time the mechanism is accurately described, but the results are limited to the fluid response, the solid wall being infinitely rigid and cannot be deformed, and consequently damaged.

Choice of numerical simulations

Each previous analysis treats separately a different aspect of the physical phenomenon, but the link between these two kinds of analysis is lacking. Moreover, experimental studies do not allow determining accurate values of physical quantities inside the droplet and the solid target. This is the reason why numerical simulations are chosen for this study. They allow calculating all physical quantities included in the model in both fluid and solid at every location and every time, depending on the discretization. Then, the physical problem deals with impact of a small area, which means that mechanical behaviour of the system is observable during a short time and on a small location.

The last reason for choosing numerical simulations is the versatility and adaptability of this method. Indeed, different physical cases can be studied just by changing the physical model as well as its parameters, like geometry or material definition.

This thesis work is organized as follows : The first chapter explains the erosion mechanism by presenting the existing studies on this damage. This will allow characterising the class of erosion due to droplet impingement and determine the strategy to simulate the corresponding damage. The second one describes the governing equations and numerical method used for the simulations themselves. The modelling choices is determined in order to simulate properly and optimally the desired phenomenon without unnecessary heavy computations. A first transient simulation is performed, and analysis of physical quantities allows understanding deeper the phenomenon. The third chapter concerns the fracture induced by mechanical fatigue (which is the erosion mechanism). General concepts are presented and the choice of a suitable fatigue criterion is investigated, this criterion allowing determining the damage intensity in simulations. Finally, a parametric analysis is performed for several droplet velocities and diameters in the last chapter, in order to determine the effect of each parameter on the damage features. Then, the effect of presence of a coating layer is investigated, as it should improve the bucket resistance.

Chapter 1

Physical phenomenon of droplet impact erosion

This part explains the physical phenomena involved inside droplet and structure responsible for the mechanism leading to the wear of metallic structures by water droplet impingement. Let us consider a spherical water droplet impacting a plane structure orthogonally. The radius of the droplet is R and it impacting velocity is V . After impact, a shock wave starts moving inside the liquid medium at the velocity c f from the contact zone and propagates also along the droplet lateral free surface. Those free surfaces act as a mirror trapping the shock. This wave follows an unobservable triple point, near the contact edge. Propagating shock front and regions under compression/tension are shown in Figure 1.3 for a perfect fluid, R = 0.1 mm V = 500 m.s -1 and ρ 0 f = 1000 kg.m -3 . The density map allows observing shock creation, propagation and interaction with the free surface. It appears clearly that free surfaces act like mirrors for waves and make them reflected perfectly. This is due to the high difference of wave velocity between the droplet and the surrounding air. After reaching the top of the droplet (region opposite to the impact area), shock wave is naturally reflected, but the superposition of relaxations waves drops the density after passing and generates tensile stress (Field [FIE 12]). Views (g) and (h) of Figure 1.3 exemplify this phenomenon. More informations about droplet impact on shock envelope droplet free surface contact edge

Contents

c L = Ì λ + 2µ ρ s c T = Ê µ ρ s (1.1) λ = νE 0 (1 + ν)(1 -2ν) µ = E 0 2(1 + ν) (1.2) 1.1.2 Contact pressure 1.1.2.

Water-hammer

During the impact of a fluid body on a solid target, the "water-hammer" pressure p wh emerges at the center of the contact area by inertial effect. Its value depends on density

p wh = ρ f c f 1 + Ȳ V Ȳ = ρ f c f ρ s c s (1.3)
Influence of Ȳ was investigated by Chen, Field and Sanada [CHE 05, FIE 89, SAN 10]. Small values of Ȳ correspond to rigid target, thus Field, Heymann, Kennedy and Li [FIE 12, HEY 69, KEN 00, LI 11] give it expression for a rigid solid body in equation (1.4), where c 0 f stands for initial sound velocity in water at room temperature :

p wh = ρ 0 f c f V (1.4)
Haller, Heymann and Li [HAL 02b, HEY 69, LI 11] approximate c f in these conditions with equation (1.5) where k f is a liquid-dependant constant, whose value equals approximatively 2 for water ([COC 97, LYO 92] for more details), v ⊥ f the particle velocity normal to the shock front, and c 0 f ≈ 1647 m.s -1 :

c f ≈ c 0 f + k f v ⊥ f (1.5)
The approximate value v ⊥ f ≈ V is taken, which gives another expression for p wh with equation (1.6). The value k f = 2 is relevant for an impact velocity range V up to 1000 m.s -1 according to Heymann [HEY 69] :

p wh = ρ 0 f c 0 f V 1 + k f V c 0 f (1.6)
The first term represents the classical water-hammer pressure derived from momentum considerations, and the second one denotes the variant nature of the shock speed (Ahmad [AHM 09]). It can be noted that the magnitude of impact pressure is independent of droplet size. However, its duration depends on droplet size and geometry according to Thomas [THO 70]. The water-hammer pressure from equation (1.6) is plotted as a function of impact velocity on Figure 1.5. The water-hammer being proportional to the square of impact velocity (equation (1.6)), its value can reach high levels for increasing velocities. For instance, V = 100 m.s -1 and 500 m.s -1 generate the corresponding values : p wh = 180 MPa and p wh = 1.3 GPa. Kennedy [KEN 00] gives the expression to evaluate radius of contact area on which the water-hammer pressure occurs R wh in equation (1.7) :

R wh = RV c f (1.7)
The first information provided by equation (1.7) is the proportional increasing dependency of R wh with droplet radius R for constant impact velocity and identical liquid. Introducing expression of c f from equation (1.5) into equation (1.7) gives the ratio bet-ween R wh and R in equation (1.8). This ratio is plotted in Figure 1.6 as function of impact velocity V .

R wh R = V c 0 f + k f V (1.8) 0
R wh /R [-] V [m.s -1 ] Figure 1
.6: Ratio between water hammer edge and droplet radius vs. impact velocity for water droplet on rigid target with c 0 f = 1647 m.s -1 .

For same droplet size, water-hammer pressure area extends while impact velocity increases, because of the inertial source of this pressure. The increase is quite linear for relatively low velocities but presents a damping effect for high velocities. Indeed, when V → ∞ we get R wh /R → 1/k f , so for water R wh → R/2. In other words, the area on which water-hammer pressure acts shall not overtakes the quarter of the droplet diameter. This shows the water-hammer pressure acts always on a narrow part of the contact zone, regardless of the intensity of impact velocity, due to the inertial aspect of water-hammer pressure.

Pressure evolution and distribution

According to Haller and Li [HAL 02b, LI 11], when a small diameter and high velocity water droplet impacts a rigid flat target, viscous effects and surface tension can be 1. Physical phenomenon of droplet impact erosion neglected. Indeed, for a droplet radius R = 0.1 mm and initial velocity V = 500 m.s -1 , Reynolds number is Re = 50 000 and Weber number is We = 350 000 [HAL 02b]. Numerical results from Haller [HAL 02b] show almost constant temperature, so convective heat transfer is not involved in the fluid motion. Therefore, the fluid behaviour is driven by inertial effects. Numerical results from Haller and Li [HAL 02b, LI 11] show the pressure distribution following the contact area across the time. According to Field, Haller, Heymann and Lesser [FIE 12, HAL 02b, HEY 69, LES 95], the maximum pressure p max occurs exactly on the edge of the contact area. Figure 1.7 illustrates the previous assumption by showing the distribution spreading along radius with a peak on the contact edge. The peak owns clearly it maximum value on intermediate time. The moment the maximum value acts is not at the start of impingement, but when jetting by lateral ejection of the fluid (Figures 1.9. The maximum pressure locates at the jetting region. These two informations are contained in equation (1.9) where R jet and T jet are respectively the location and the time of jetting :

p max = p(x = R jet ,t = T jet ) (1.9)
The jetting phenomenon is detailed in the next section (1.1.3). Unfortunately analytical expression for maximum pressure does not exist. Numerical results of Haller and Kennedy [HAL 02b, KEN 00] give respectively p max 2p wh and p max 3p wh . In his theory, Cook [COO 28] shows that pressure generated at liquid-solid impact is sufficient to cause erosion in steam turbine blades. Field [FIE 12] proposes an expression on the time after which pressure realease starts T p . Release waves reach the central axis (x = 0) and terminate the high-pressure stage after the time T hp given by Bowden [BOW 64]. These two quantities are written in equation (1.10) :

T p = RV 2c 2 f T hp = 3RV 2c 2 f (1.10)
The high pressure duration T hp is then devlopped with c f = c 0 f + k f V from equation (1.5), which gives equation (1.11) :

T hp = 3RV 2(c 0 f + k f V ) 2 (1.11)
High pressure duration increases proportionally to droplet radius R, its dependency to impact velocity V is different. The time starts to increase with V , then deacreases. Its maximum value is found for ∂T hp /∂V = 0, and corresponds to V = c 0 f /k f . Note that it only depends on fluid constants except the droplet size. For water (see section 1.1.2.1 for k and c 0 f values), this velocity corresponds to V ≈ 800 m.s -1 . T hp is graphically represented in Figure 1.8 as a function of R and V . This map shows that T hp < 120 ns, even for the largest droplets at R = 1 mm. This very short time proves the intensity of the shock involved with droplet impact.

Once incompressible stream line flow is established, the pressure on the central axis falls to the much lower Bernoulli stagnation pressure p sta given in equation (1.12) ([FIE 12]) :

p sta = 1 2 ρ 0 f V 2 (1.12)

Jetting

When the contact edge becomes subsonic, shock wave overtakes it. At this exact moment, compression with solid leads to jetting by lateral ejection of the fluid (Figure 1.9). Figure 1.10 shows lateral velocity and density distribution of the fluid on contact area at the exact time of jetting. The jetting position R jet is found where density roughly decreases, while velocity presents a clean extremum. For x < R jet , the liquid is compressed and the density is quite constant and superior to the initial value ρ f > ρ 0 f , and diminishes sharply at x = R jet to finally vanish, because of the phase change, the zone x > R jet being localized outside of the liquid core. The lateral velocity is zero on the center of the droplet x = 0, because of the symmetric and continuous behaviour of the droplet, and increases slowly then promptly around x = R jet , then decreases, but holds a hump, representing air expelled from the gap between droplet and substrate. The velocity of the jet can be far higher than the impact velocity V and even the ambient sound velocity c 0 f . Indeed, Figure 1.10 shows a jet velocity higher than 2000 m.s -1 , whereas the sound velocity is here c 0 f = 1647 m.s -1 . The jetting time T jet is defined as the time when the liquid medium breaks (spalls) through the droplet free surface at the contact line. From a theoretical consideration, this can be expected to occur when the contact line velocity becomes equal to the shock velocity at the contact line. Haller [HAL 02b] suggests the time when jets form with equation (1.13), where ĉ f stands for compression wave velocity inside the droplet when jetting :

T jet = RV 2 ĉ2 f (1.

Cavitation

Far-side cavitation

When the main shock wave is reflected opposite to the impact point, the spherical nature of the free surface forces the corresponding reflected wave (labelled 'R' in Figure 1.12) to converge towards an unique point or "focus" (labelled 'F' in Figure 1.12). Cavitation can be produced in this region according to Haller, Obreschkow, Sanada and Xiong [HAL 02b, OBR 11, SAN 08, XIO 10]. When droplet is spread out later, this cavitation zone can be close enough to the solid to cause damage by bubble collapsing. This kind of cavitation is only active in some particular cases such as aircraft and missiles eroded by rain [FYA 66, ZAH 81] and Pelton turbine buckets eroded by highspeed (up to 200 m.s -1 ) droplets and jets [PER 07]. Experimental results of Obreschkow [OBR 11] show the absence of erosion due to far-side cavitation for V = 110 m.s -1 and R = 5 mm, because focal region appears too far from contact area. The roughening at the interface is also caused by cavitation.

Erosion mechanism : experimental approach 1.2.1 Test rigs

One of the most used experimental equipments for the study of droplet impact erosion is the Pulsating Jet Erosion Test Rig (PJET), based on high velocity jet emission as developed by Cavendish Laboratory [FIE 99]. At first, focused water jet is generated normal to the sample by a high pressure pump through a nozzle. It is then divided into coherent segments by a rotating disc with orifices [LUI 13]. Figure 1.13 illustrates its operating principles, and more details are given by Lammel and Tobin [LAM 12, TOB 11], who use the PJET in their inversigations.

The second one is the Rotating Arm Appartus. This rig is especially used to reproduce rain impact with droplets. In this device, one or more specimens are attached to the periphery of a rotating disk or arm, and their circular path passes through one or more liquid jets or sprays, causing discrete impacts between the specimen and the droplets or the cylindrical surface of the jets [ AST 10 (i) The first phase is an initial region, called "incubation period". During this time no significant loss in mass is observed, but the surface condition changes and becomes more rough.

(ii) A region in which the rate of erosion increases steadily (iii) A region in which the rate of erosion is at a maximum and a sensibly constant (iv) Finally, erosion rate diminishes, possibly again becomes constant, or zero in some cases.

These steps are illustrated in Figure 1.16. Ma [MA 15] splits the mechanism into 5 parts. The only difference concerns the last step, which is divided into a deceleration phase, then a termination one. Baker [BAK 66] proposes analytical expression of the maximum erosion rate ξ max in equation (1.14) which corresponds to the slope of erosion curve on phase (iii) shown in Figure 1.16 :

ξ max = ∂M e ∂q e max = A (V sin θ -V c ) n 1 csx θ (1.14)
where M e and q e stand respectively for the mass loss and mass of water impacting per unit surface. The angle of incidence is described by θ, thus, the term V sin θ represents the normal velocity. A, V c and n 1 are constants characteristic of the material. V c may be interpreted physically as the value of the impact velocity below which erosion will not occur and in this sense it is analogous to the fatigue limit in fatigue tests. Therefore, equation (1.14) is interpreted in terms of a fatigue mechanism damage. This relation is given for a constant diameter. The corresponding value of n 1 is 4 -5 for ductile materials and 6 -9 for brittle materials (Ahmad [AHM 09]). The equation shows a maximal value for θ = 90 • , i.e. for a perpendicular impact. Staniša [STA 92] studied the influence of impact angle on the steam turbine blades erosion. According to Hancox [HAN 66], among other parameters, impact velocity V , impacting droplet size, or diameter φ = 2R, impacting water quantity q e and impacting angle θ are the main parameters which can significantly change the erosion pattern in any testing environment. Oka [OKA 07] assures the importance of impact frequency for damage erosion and proposes an enhanced expression for the damage definition.

Another model was proposed by [ SAN 88] and used by [FUJ 12] to evaluate the volume loss inside nuclear power plant pipes and takes into account the hardness of the target and considers the cumulative plastic strain.

Influence of impact angle

The influence of θ on erosion rate was studied by Ahmad and Hattori [AHM 09, HAT 13]. An example of result is presented on the Figure 1.17, which corroborates the equation (1.14). This graphic results from experimental investigations on X20Cr13, material used as a blade steel, and gives the volume loss for different angles for the same exposure duration, which time evolution has the same shape as erosion rate of phase (iii) ξ max (θ) given in equation (1.14). Finally, for a normal droplet impact, equation (1.14) can be simplified and put in the form of equation (1.15) : M e ∝ φ n 2 (1.16)

ξ max (θ = 90 • ) ∝ V n 1 (1.
Figure 1.18 shows the evolution of volume loss as a function of number of impacting droplets, relative to their size, in unit per cm 2 . One can notice that erosion threshold starts later for smaller droplets. In other words, incubation phase defined in section 1.2.2.1, depends on droplet size. Moreover, the slope of the curve is higher for large droplets, thus the gap between small and large droplets widens with time.

Influence of impact velocity

Influence of impact velocity V was studied by Gerdes, Oka and Mahdipoor [GER 95, OKA 07, MAH 15] for steam turbine blades and by Li [LI 12] for nuclear power plant pipes. Gerdes [GER 95] compared resistance of several materials and effects of laser nitrining. The relation between incubation period and droplet velocity is similar to an S-N curve for fatigue. The relation between the impingement velocity and erosion damage per unit mass of water droplets indicated a high exponent number of about 6 in the region of low velocity because the water droplet intensity against the material strength rapidly decreases with the decrease in droplet velocity. However, the velocity exponent of 2 was recognized in the high velocity region. This indicates that the erosion damage is basically related to the energy of impinging droplets in the high velocity region. This is similar to erosion by solid particle impact [OKA 07]. Figure 1.19 shows the cumulative volume loss as a function of number of impacts for several impact velocities. Like droplet size, velocity has a huge influence on incubation period, where erosion starts earlier for high velocity, and the slope of the curve is higher too for high velocities. 2. These regions undergo plastic strain, which produces residual compression stresses. This underground plasticity is responsible for surface shape change. Slip bands can be detected on surface (Figure 1.25).

3. High strain cycles are generated and cracks form around grains boundary (Figure 1.23).

4. These cracks grow, then link together, and form a underground cracks network.

5. After high number of deformation cycles, grain tilting is observed. Moreover, intergranular cracks born, producing microvoids on the surface (Figure 1.24).

6. When intergranular cracks networks are developed enough, gains are ejected and reveal triple junctions, unleashing pits on the surface (Figure 1.25). Striations locate inside these pits, proving the cyclic nature of the damage, i.e. fatigue cracking. Kong [KON 10] proposes a more suitable description for grain-sized droplets, and split them into initiation and propagation phases, like macrospocopic steps listed in previous section (1.2.2.1). These steps are listed below :

Physical phenomenon of droplet impact erosion

→ Initiation and development

Step 1 : Plastic deformation is generated and induces crack initiation. For brittle material, crack lines are propagated outwards from the initial contact area.

Step 2 : Propagation of compression waves produces tensile stresses, which act across the circumference of the deformation zone in high velocity impact Step 3 : Microscopic pits are produced by cracks intersection, and are step-like opening. Then surface layer seems to be torn away. Microcracks at deformation zone join together, which generates small zone of removed material. This kind of damage is observed for cavitation too [HAO 08], but erosion shape is smaller and more oval.

→ Growth and propagation

Step 4.a : Intergranular cracks appear and are observable in surface, thus fracture surface has the appearance of grain structure, with secondary cracks [BLA 85]. Intergranular cracks tend to open gap between grains (Figure 1.27).

Step 4.b : Triple junctions emerge due to structural imperfections. Indeed, crystal lattices are randomly oriented, and structure contains dislocations, voids, etc... These imperfections are responsible for weakening of grain boundaries (Figure 1.28).

Step 4.c : Intragranular cracks appear in addition to intergranular ones. Grains are fractured in two modes : translamellar and interlamellar, depending on the lamellar orientation (Figures 1.27.a and 1.29).

Coating

All mechanisms described in the previous sections are responsible for erosion of Pelton runners, in addition to others mechanisms not considered in this work. The vulnerable zones of a Pelton bucket are marked in Figure 1.30. In order to reduce the wear of buckets subject to liquid impact, manufacturers improve their resistance with surface treatments. One of the most known in hydraulic machinery 1. Physical phenomenon of droplet impact erosion -High velocity oxy-fuel coating spraying (HVOF)

-High velocity air fuel (HVAF)

-Warm spraying -Cold spraying Among all these technologies, most of hydraulic turbines manufacturers turn to High velocity oxy-fuel coating spraying , or HVOF. An example of Pelton runner being coated by HVOF process in shown in Figure 1.31. This process is a relatively recent addition to the family of thermal spray processes. During the 1980s, a class of thermal spray processes called high velocity oxy-fuel spraying was developed. Figure 1.32 describes the general principle of HVOF process. A mixture of gaseous or liquid fuel and oxygen is fed into a combustion chamber, where they are ignited and combusted continuously. The resultant hot gas at a pressure close to 1 MPa emanates through a converging-diverging nozzle and travels through a straight section. The fuels can be gases (hydrogen, methane, propane, propylene, acetylene, natural gas, etc.) or liquids (kerosene, etc.). The jet velocity at the exit of the barrel (> 1000 m.s -1 ) exceeds the speed of sound. A powder feed stock is injected into the gas stream, which accelerates the powder up to 800 m.s -1 . A illustration of a couple of coating powders materials is proposed in Figure 1.33. The stream of hot gas and powder is directed towards the surface to be coated. The powder partially melts in the stream, and deposits upon the substrate (Figure 1

Physical phenomenon of droplet impact erosion

Conclusion on the mechanism of droplet impact erosion

This chapter described the mechanisms responsible for erosion of metallic structure, as Pelton buckets by water droplet impingement. The first part deals with the fluid-structure interaction itself between the spherical water droplet and the plane solid target. Concerning the fluid, at the beginning of impact, a shock wave travels from the contact zone to the opposite side. The free surface acts like a mirror and generates a reflected wave in the fluid. After reaching the top of droplet, the superposition of relaxation waves drops the pressure. The compression waves are present inside the solid in the same way, starting from the impact point, but are followed by shear waves too, and to a lesser extend by Rayleigh waves on the surface. The contact pressure can be firstly described by the waterhammer pressure, whose value can be analytically found, and denotes the pressure on the center of the contact at start of impingement. The maximum contact pressure, higher than water-hammer pressure, locates on the edge of contact area, and is found when microjets form. After that, its value diminishes to the stagnation pressure. These microjets are generated by the compression effects of the shock wave travelling inside the fluid. When the expanding contact area becomes subsonic, shock front overtakes it and compression effects eject the liquid in contact with solid. The velocity of microjets is far higher than the ambient sound speed of liquid. The reflected waves inside the droplet force local density to vanish, and this phenomenon could induce cavitation by phase change, but this behaviour will not be included to forthcoming simulations. For small droplet at high velocity, fluid behaviour is driven by inertial effects, thus viscosity and surface tension could be neglected for a relevant simulation.

The second part concerns the erosion mechanism itself, by experimental approach found in literature. The main test rigs for droplet impact erosion are presented, such as the Pulsating Jet Test Rig (PJET) or the Rotating Arm Apparatus, which is often used for Conclusion on the mechanism of droplet impact erosion rain impact erosion. Then, the macroscopic models and experimental studies on erosion are described. These models and studies give a general description of different phases observed : an incubation phase, followed by a linear erosion, and finally the last one where erosion rate decreases. Experimental investigations from literature analyse the influence of parameters like droplet velocity, its size or impact angle. It has been found that high velocity and size with normal impact increase damage. Next, the mesoscopic observation of the damage is presented through experimental works from literature. The mechanism is described by different steps. First, the first impacts erode grain boundaries and generate pits between grains. Then, microckacks appear at the bottom of these pits. After a larger number of impacts, material is removed from surface when cracks join together. The erosion mechanism is finally fatigue cracking according to both macroscopic and mesocopic observations. This last information will be basic mechanism for the development process of the erosion tool in Chapter 3.

The last part presents the surface treatment like coating, used to increase erosion resistance of hydraulic machines. The High velocity oxy-fuel coating spraying, or HVOF, which is widely used for Pelton buckets is exposed.

Chapter 2 Numerical methods for droplet impact

This second chapter gives an overview of physical and numerical equations acting in the considered system. Then, the numerical convergence of models is checked and one case of coupled computation is presented. This section presents the analytical equations of continuum mechanics driving the fluid and solid behaviour. First the general continuum equations are described, and then are translated in terms of solid and fluid, with the constitutive laws and state equations. Finally, the coupling conditions between the fluid (water droplet) and the solid (bucket) sub-domains are specified.

Different configurations in continuum mechanics

Let us consider a body by its material domain Ω = Ω(t), bounded by a closed surface ∂Ω = ∂Ω(t). Its initial state is defined with Ω 0 = Ω(t = 0) and ∂Ω 0 = ∂Ω(t = 0). The reference (or undeformed) configuration is referred to C 0 and given by (Ω 0 , ∂Ω 0 ) and the current one C (t) by (Ω, ∂Ω). The current position in C 0 is denoted X = X(t), and is located by x = x(t) in C (t). The variable X is called material coordinates or Lagrangian coordinates, and x is called spatial coordinates or Eulerian coordinates. The transformation between the two configurations is noted Φ and given in equation (2.1) :

x = Φ(X,t) X = Φ -1 (x,t) (2.1)
If a function is expressed as F(X,t), it depends on the Lagrangian coordinates X and the time t : this is called the Lagrangian description. The function is defined at the initial coordinates which are not time-dependent. Whereas, if the function is expressed as f (x,t), it will depend on the spatial or Eulerian coordinates x and time t. This function is then defined at the current coordinates which are time-dependant, and it is called the Eulerian description. More informations about these descriptions are provided by [LET 09, LI 13]. To avoid the drawbacks of each of the two classical descriptions, a more generalized description method is needed, which is called Arbitrary Lagrangian-Eulerian (ALE) approach. In the ALE description of motion, neither the materiel configuration, nor the spatial configuration is taken as the reference. The physical variables and equations are expressed in a third domain which is named as the referential configuration where the coordinates are denoted by χ. With another application Ψ, this reference coordinate χ is related with the material coordinates X, and with the spatial corrdinate x by Ò Φ. The motion of the computational mesh is described by equation (2.2) :

χ = Ψ(X,t) x = Ò Φ(χ,t) (2.2)

Conservations laws in continuum mechanics

The continuum mechanics assures the conservations of three important quantities along time inside a domain called Ω, bounded by ∂Ω. These quantities are : mass, momentum and energy and their conservation laws are obtained by the Reynolds transport theorem (RTT) given in equation (2.3), where f (x,t) is the desired quantity, n the normal to ∂Ω and v(x,t) the velocity :

∂ ∂t X Ω f (x,t) dx = Ω ∂ f (x,t) ∂t x dx + ∂Ω f (x,t)v(x,t) • n dx (2.3) Conservation of mass Taking f (x,t) = ρ(x,t) in equation (2.
3), which leads to equation (2.5) with Gauss's theorem :

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

∂ρ ∂t x + ∇ x • (ρv) = 0 (2.5) Conservation of momentum Taking f = ρv in equation (2.
3), which leads to equation (2.6) :

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

Ω ∂ρv ∂t x dΩ + ∂Ω (ρv ⊗ v) • n dΓ = Ω ρb dΩ + ∂Ω σ • n dΓ (2.7)
where b denotes the body force vector and σ the Cauchy stress tensor. After using the Gauss's theorem, the conservation of momentum is given in equation (2.8) :

∂ρv ∂t x + ∇ x • (ρv ⊗ v) = ρb + ∇ x • σ (2.8)

Conservation of energy

Taking f = ρe 0 , where e 0 = e + 1 2 ρv • v is the total energy in equation (2.3), which leads to equation (2.9) :

∂ ∂t X Ω ρe + 1 2 ρv • v dΩ = Ω ∂ ∂t x ρe + 1 2 ρv • v dΩ + ∂Ω ρe + 1 2 ρv • v v • n dΓ (2.9)
where e denotes the internal energy per unit volume. The material time derivation of the total energy over the material domain Ω equals the power of the volume and surface energy sources. Hence, the balance of the energy can be expressed by equation (2.10) :

Ω ∂ ∂t x ρe + 1 2 ρv • v dΩ + ∂Ω ρe + 1 2 ρv • v v • n dΓ = Ω ρb • v dΩ + ∂Ω t n • v dΓ + Ω ρs dΩ - ∂Ω q • n dΓ (2.10)
where the stress vector t n = σ • n denotes the contact force exerted on the boundary of the material domain, ρs the heat source per unit volume, and q the heat flux per unit area. After using the Gauss's theorem, energy balance takes the form of equation (2.11) :

∂ρe 0 ∂t x + ∇ x • (ρe 0 v) = σ : L -∇ x • q (2.11)
where L = ∇ x v denotes the velocity gradient tensor, defined as L i j = ∂v i ∂x j . As mentioned in section 1.1.2.2, the interaction between the solid and fluid consists only in the exchange of linear momentum (inertial effect), hence the thermal phenomena are neglected, i.e, there is no heat transfer at the interface, or heat creation inside the material bodies. As a consequence, in the present work, convective effects shall be ignored and equation of energy unused.

The Solid sub-domain 2.1.3.1 Conservation equations

Eulerian description

The conservations equations (2.5) and (2.8) are adapted to solid version in equation (2.12) :

         dρ s dt + ρ s ∇ x • v s = 0 ρ s dv s dt = ∇ x • σ s + ρ s b (2.12)
This system of governing equations in Eulerian description is usually used with the Finite Element Method (FEM) in Update Lagrangian Fomulation (ULF). One important advantage of this last formulation is that the stress and strain are defined in the spatial domain, which means that they both have the intuitive physical senses. Hence the constitutive models are easily to be implied in the numerical simulations. As presented previously, a physical variable can also be expressed with the Lagrangian description approach, i.e. with the material coordinate X.

Lagrangian description

Here the governing equations in Lagrangian description are given directly in equation (2.13) :

       ρ s (X,t)J(X,t) = ρ 0 s (X) ρ 0 s ∂v s (X,t) ∂t X = ∇ X • P + ρ 0 s (X)b (2.13)
where ρ 0 s (X) denotes the initial solid density which depends on only the material coordinates X, J(X,t) is the determinant of the Jacobian matrix F = ∇ X x of the motion Φ(X,t). The matrix P = Jσ • F -is the first Piola-Kirchhoff stress tensor and is expressed in the initial configuration Ω 0 . The nominal stress tensor N = P is sometimes used instead. More details can be found in [BEL 13].

Constitutive relation

The stress tensor σ s has to be expressed by a constitutive equation as a function of the strain ε s . The material will be small strains elasto-plastic, thus the strain can be splitted into an elastic part ε e s and a plastic one ε p s in equation (2.14) :

ε s = ε e s + ε p s (2.14)
The strain is a function of the displacement u s whose material derivative is the velocity v s = ∂u s ∂t X in equation :

ε s = 1 2 ∇ x u s + ∇ x u s + ∇ x u s • ∇ x u s (2.15)
The elastic part of strain is linked to stress by the the stiffness tensor C in equation (2.16) in Eulerian description :

σ s = C : ε e s (2.16)
This relation is also called Hooke's law. In case of a linear elastic and isotropic material, the relation (2.16) simplifies and only depends on Lamé coefficients (λ, µ) defined in section 1.1.1.2 (Young's modulus E and Poisson's ratio ν are sometimes used instead of λ and µ). The new form is given in equation (2.17) :

σ s = λtr (ε e s ) 1 + 2µε e s (2.17)
where 1 denotes the unit tensor. The corresponding Lagrangian description of relation (2.16) is S = C : E, where S = JF -1 • σ • F -denotes the second Piola-Kirchhoff stress tensor and E = 1 2 (F • F -1) stands for the Green-Lagrange strain tensor. The use of S is justified by its property of symmetry, unlike P.

The general strain (including pasticity) is linked to the stress by the equation (2.18) :

   dσ s = H ep (σ s , X) : dε s F Y (σ s , X) = σ VM -(σ 0 Y + X) (2.18)
where H ep (σ s , X) is the elastoplastic tangent modulus, F Y (σ s , X) stands for the yield function, σ 0 Y denotes the yield strength, and X is the isotropic hardening variable, which depends on the cumulative plastic strain εp . The Von Mises stress σ VM is chosen as equivalent stress for plasticity criterion. These two last quantities are described below. More informations about constitutive models are found in [AUN 16,LEM 94].

All these equations constitute the governing equations for the solid sub-domain, which are to be discretized in space and integrated in time.

Quantities derived from governing equations

The Cauchy stress σ s is a symmetric second-order tensor, hence it can be splitted into a spherical part σ sp = σ H 1, where σ H is called hydrostatic stress and a deviatoric one s, the deviatoric stress :

σ s = s + σ H 1 (2.19)
where

σ H = 1 3 tr (σ s ) s = σ s -σ H 1 (2.20)
The spherical part is a diagonal matrix with a constant value on its diagonal, whereas the deviatoric one is a zero-trace tensor. In the same way, the strain tensor ε s can be splitted into a spherical part, called hydrostatic strain ε H , and a deviatoric one e defined in equation (2.21) :

ε H = 1 3 tr (ε s ) e = ε s -ε H 1 (2.21)
The total strain energy per unit volume is defined in equation (2.22) :

U = 1 2 σ s : ε s (2.22)
Using the splitting defined in equations (2.20) and (2.21), the strain energy is itself splitted into a part relative to change of volume U v and a part relative to distortion or change of shape U d such as U = U v +U d in equation (2.23) :

U v = 3 2 σ H ε H U d = 1 2 s : e (2.

23)

The deviatoric stress, which mesures the difference between the stress vector and mean normal stress, is an important quantity for fracture mechanics and the study of plasticity. The Von Mises stress σ VM is proportional to the euclidean norm of s and is defined in equation (2.24) :

σ VM = Ê 3 2 s : s (2.24)
Being a norm the Von Mises stress is a positive quantity. Let us define the signed Von Mises stress σ sign VM by adding the sign of the trace of stress in equation (2.25) :

σ sign VM = sign (tr (σ s )) σ VM = σ H |σ H | σ VM (2.25)
The Von Mises equivalent plastic strain rate εp is defined as an euclidean norm of the plastic strain rate tensor εp = ∂ε p s ∂t in equation (2.26) :

εp = Ê 2 3 εp : εp (2.26)
The cumulative or equivalent plastic strain εp at the time t, is then the time-integral of the equivalent plastic strain rate in equation (2.27) :

εp (t) = t 0 εp (u) du (2.27)
The volumetric dilation Θ s is linked to the hydrostatic strain, and for infinitesimal strains is the divergence of the displacement :

Θ s = 3ε H ≈ ε 1 ∇ x • u s (2.28)
For linear elastic and isotropic material, hydrostatic stress σ H is linked to hydrostatic strain ε H by the bulk modulus K s . By using equation (2.28) one gets a linear relationship between hydrostatic stress and dilation :

σ H = K s ε H = K 3 Θ s K s = 3λ + 2µ (2.29)
The stress triaxiality ζ is defined as the ratio between hydrostatic and Von Mises stresses :

ζ = σ H σ VM (2.30)
Its magnitude and sign give informations about the nature of load occurring inside the structure at given time and location. Some particular values are listed below :

• ζ = 0 : pure shear • ζ = 0.33 : uniaxial traction • ζ = 0.66 : equi-bitraction • ζ = -0.33 : uniaxial compression • ζ = ±∞ : hydrostatic pressure 2.1.3.

Invariants and principal stresses

The Cauchy stress tensor σ ≡ σ i j being symmetric, it can be diagonalized. There are three eigenvectors, called principal vector or axis e I , e II , e III associated to three eigenvalues, called principal stresses σ I , σ II and σ III , respectively. The stress tensor having any form in reference basis B r = (e 1 , e 2 , e 3 ), becomes a diagonal matrix when expressed in its principal basis B p = (e I , e II , e III ) :

σ = σ 11 σ 12 σ 13 σ 22 σ 23 sym σ 33 B r = σ I 0 0 0 σ II 0 0 0 σ III B p (2.31)
Expression (2.31) shows that every plane formed by principal axis supports only traction-compression, without shear. Indeed, the stress vector acting on a plane normal to any principal direction e L , with L = I, ..., III is collinear to e L itself and its intensity equal to the corresponding principal stress σ L :

t e L = σ • e L = σ L e L (2.32)
Moreover, the highest principal stress, in absolute value, is the maximum physical stress inside the material at the current time and point, and is often used as plasticity/fracture criterion, like Rankine criterion.

The principal stresses are the solutions of the characteristic equation : det (σλ σ 1) = 0, where λ σ = σ I...III denotes the set of eigenvalue. This allows defining the characteristic polynomial P σ :

P σ = -λ 3 σ + I 1 λ 2 σ -I 2 λ σ + I 3 (2.33)
where I 1 , I 2 , I 3 are the invariants of σ. They do not depend on the basis in which the stresses are expressed and are given by :

               I 1 = tr (σ) I 2 = 1 2 tr (σ) 2 I 3 = 1 3 tr (σ) 3 = det (σ) (2.34)
The same procedure can be applied to the deviatoric stress s, whose characteristic polynomial P s is, with λ s as set of eigenvalues : 2. Numerical methods for droplet impact

P s = -λ 3 s + J 1 λ 2 s -J 2 λ s + J 3 (2.35)
where J 1 , J 2 , J 3 are the invariants of s. They can be expressed as functions of stresses invariants :

               J 1 = tr (s) = 0 J 2 = 1 2 tr (s) 2 = 1 3 I 2 1 -I 2 J 3 = 1 3 tr (s) 3 = det (s) = 2 27 I 3 1 - 1 3 I 1 I 2 + I 3 (2.36)
Note that s has the same principal axis as σ, because expressing σ as a function of s in its characteristic equation P σ = 0 leads to : det (σλ s 1) = 0, with λ s = λ σσ H , i.e. :

       s I = σ I -σ H s II = σ II -σ H s III = σ III -σ H (2.37)
We can observe that hydrostatic and Von Mises stresses can be expressed as functions of these invariants :

     σ H = I I 3 σ VM = √ 3J 2 (2.38)
Equation (2.38) shows that hydrostatic and Von Mises stresses do not depend on the basis, which makes them very convenient for stress analysis.

The principal axis allow building the octahedral planes. They are equally inclined to the directions of the three principal stresses and form an octahedron. The specificity of this octahedron is that each of its faces has the same normal and tangential stresses. The stresses acting on these planes are called normal octahedral stress σ oct and tangential or shear octahedral stress τ oct and can be simply expressed as functions of the previous stresses :

     σ oct = σ H τ oct = √ 2 3 σ VM (2.39)
Using the principal stresses, we build the principal stresses ratios r i, j defined as : Ω is moving with an arbitrary velocity v. The stress tensor inside fluid σ f can be simplified here. Indeed, the Reynolds number applied to the present droplet impact case can be set into the form of equation (2.41) proposed by [HAL 02b] :

r 1,2 = σ I σ II r 2,3 = σ II σ III r 3,1 = σ III σ I (2.40)
Re = ρ 0 f V R η (2.41)
where η stands for the dynamic viscosity of the fluid. For the present water droplet, ρ 0 f = 1000 kg.m -3 , V = 100 m.s -1 , R = 0.5 mm and η = 9 • 10 -4 Pa.s, we get a Reynolds number of Re = 55 000, which shows the predominance of inertial effects over viscous ones. Hence, σ f does not contain viscous term, but only the hydrostatic one and becomes a spherical tensor in equation (2.42) :

σ f = -p f 1 (2.42)
The fluid is then perfect, or inviscid, and Euler equations shall be used as simplification of Navier-Stokes ones to compute its mechanical behaviour.

The developments lead to the final system of equations (2.43) (see [LI 13] for more details) :

             ∂ ∂t χ Ω ρ f dΩ + Ω ∇ x • ρ f (v f -v) dΩ = 0 ∂ ∂t χ Ω ρ f v f dΩ + Ω ∇ x • ρ f v f ⊗ (v f -v) + p f 1 dΩ = Ω ρ f b dΩ (2.43)
In addition, the fluid flows are considered isothermal-barotropic, which means that the fluid temperature is constant and the fluid pressure is only a function of fluid density (and vice-versa). The modified Tait's equation of state proposed by [MAC 66] is chosen for water and is given in equation (2.44) :

p f = B   ρ f ρ ref f γ f -1   + p ref f (2.44) where B = ρ ref f (c ref f ) 2 /γ f . The following reference quantities are ρ ref f = ρ 0 f , c ref f = c 0 f and p ref f = 0,
and the polytropic coefficient γ f ≈ 7 is taken. Finally, the equation of state becomes :

2. Numerical methods for droplet impact

p f = ρ 0 f (c 0 f ) 2 γ f   ρ f ρ 0 f γ f -1   (2.45)
Equations (2.43) and (2.45) constitute the governing equations for the fluid subdomain, which are to be discretized in space and integrated in time.

The last important feature eventually involved in the liquid behaviour is the surface tension. The Weber number for the present application is given by [HAL 02b] :

We = ρ 0 f V 2 R ς (2.46)
where ς denotes the surface tension coefficient between the two regions : water and air in our case. For the present water droplet, ρ 0 f = 1000 kg.m -3 , V = 100 m.s -1 , R = 0.5 mm and ς = 72 N.mm -1 , we get a Weber number equal to We = 70 000, which shows the predominance of inertial effects over surface tension. Therefore, surface tension shall not be considered in the calculation of the fluid response.

Coupling conditions

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

For a viscous fluid flow interacting with a deformable structure, the continuity conditions at the interface are usually given as equation (2.47) :

∀x ∈ Γ I v f (x,t) = v s (x,t) σ f (x,t) • n f (x,t) = -σ s (x,t) • n s (x,t) (2.47)
where Γ I denotes the fluid-structure interface, n f (x,t) and n s (x,t) are the normal vectors pointing-out of each sub-domain at the interface, with n f (x,t) = -n s (x,t) as shown in Figure 2.1. One can not that the second part concerns the stress vector and is equivalent to

t f ,n f (x,t) = -t s,n s (x,t).
For an inviscid fluid, the continuities of velocity and force are not enforced in the tangential direction, hence the slipping condition or non-penetrating condition is used at the interface in equation (2.48), which imposes only the continuities in the normal direction :

∀x ∈ Γ I    n f (x,t) • v f (x,t) = n s (x,t) • v s (x,t) n f (x,t) • σ f (x,t) • n f (x,t) = n s (x,t) • [σ s (x,t) • n s (x,t)] (2.48)
Because in this work the fluid is considered to be inviscid, we will apply the interface condition of equation ( 2 presented in the previous section, the stress section can be written as σ f = -p f 1. Because of the slipping condition, the condition applied to the stress can be written otherwise with the normal stress at the interface σ s,n s (x,t) = n s (x,t) • t s,n s (x,t), which for the sake of brevity, will be noted p s = -σ s,n s , as the pressure applied on the solid at the interface. Hence, the condition at the interface Γ I becomes equation (2.49), and only deals with normal velocity and stress, or pressure : Notice that both the two formulations are with Lagrangian meshes, which means that the nodes and elements move with the material, and the boundaries and interfaces remain coincident with element edges. Furthermore, the constitutive equations are always evaluated at the same material points, which is advantageous for history-dependent materials. More concretely, the time derivatives in the governing equations (2.12) or (2.13) are both carried out holding the material coordinates constant. The differences between the two formulations are that in ULF the spatial derivative is with respect to the spatial coordinates and the corresponding weak form is obtained by integration over the current configuration of the solid body, whereas in TLF the spatial derivative is with respect to the material coordinates and its weak form is given by integration over the initial configuration. Because the ULF and TLF are mathematically equivalent, only the weak form for the ULF is presented in this section. e Ω e . Two different types of boundary conditions are imposed :

n f • v f = n s • v s p f = p s (2.
Ω s Γ t Γ v Ω e
• Velocity boundary condition v(X,t) = v b (X,t) on Γ v (2.50) • Traction boundary condition σ(X,t) • n b (X,t) = t b (X,t) on Γ t (2.51)
where v(X,t) denotes the velocity vector expressed with the material coordinates X, v b is the imposed velocity at the boundary Γ v . σ(X,t) represents the stress tensor, n b the normal vector of the boundary Γ t and t b (X,t) the traction or force exerted on Γ t . The relation between Γ v and Γ t is described as (2.52), where Γ s denotes the boundary of the whole solid domain Ω s :

Γ v ∪ Γ t = Γ s Γ v ∩ Γ t = ∅ (2.52)
The weak form is obtained by multiplying the momentum equation by a test function and then integrating over the current domain, which gives equation (2.53) :

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

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

v(X,t) ∈ U U = {v(X,t) | v ∈ C 0 (X), v = v b on Γ v } (2.
Ω s δL : σ dΩ - Ω s δv • ρb dΩ - Γ t δv • t b dΓ + Ω s δv • ρv dΩ = 0 (2.56)
where δL = ∇ x δv is the spatial gradient tensor of the virtual velocity. This tensor can be decomposed into symmetric and skew-symmetric parts in equation (2.57) :

δL = δD + δW (2.57) with δD = 1 2 (δL + δL ) δW = 1 2 (δL -δL ) (2.58)
where δD is called the virtuel strain rate, which is the symmetric part of δL, and the spin δW, called the virtual rotation rate is the skew-symmetric part of δL. Because the Cauchy stress tensor σ is symmetric, i.e. σ = σ , then we have δL : σ = δD : σ, since δW : σ = 0. The procedure for obtaining the weak form is also called the Virtual power principle, in which the the test function δv is called the virtual velocity, the virtual power δP tot of the whole system is equal to zero. In fact, in equation (2.56) each part possesses a physical sense :

• Virtual internal power δP int δP int = Ω s δD : σ dΩ (2.
δP tot = δP int -δP ext + δP kin = 0 ∀δv ∈ U 0 (2.62)
This weak form for Updated Lagrangian Formulation will be discretized in next section.

Finite Element discretization

In the Finite Element Method (FEM), the motion x(X,t) is approximated by equation (2.63) [BEL 13, ZIE 77] :

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

x(X,t) = N i (X)x i (t) (2.64)
Using the same approximation methodology to the displacement field, one obtains :

u(X,t) = N i (X)u i (t) (2.65)
where u i (t) = x i (t) -X i with X i being the initial position of the node i. The velocity is obtained by taking the material time derivative of the displacement :

v(X,t) = ∂u(X,t) ∂t X = ∂[N i (X)u i (t)] ∂t X = N i (X) ui (t) (2.66)
Similarly, the acceleration field is approximated by :

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

δv(X) = N i (X)δv i (2.68)
with which δL can be written as

δL = ∂δv ∂x = δv i ⊗ B i = δv i B i (2.69)
where B i = ∇ x N i (X) stands for the gradient of shape functions.

Remark : In this work, a variable or an equation may be expressed either in tensor form or in matrix form, which are essentially equivalent. When the operation is carried out with an operator like "•" or "⊗" etc., it means that we are using the tensor form. Whereas, if in the equation it has the transpose symbol " ", it is in the matrix form.

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

                   δP int = δv i Ω s σB i dΩ = δv i f int i δP ext = δv i Ω s N i ρb dΩ + Γ t N i t b dΓ = δv i f ext i δP kin = δv i Ω s N i ρN j dΩ v j = δv i M i j v j (2.70)
where f int i and f ext i are the internal and external nodal forces for the node i, respectively, and M i j is one of the elements in the so called mass matrix M :

                 f int i = Ω s σB i dΩ f ext i = Ω s N i ρb dΩ + Γ t N i t b dΓ M i j = Ω s N i ρN j dΩ (2.71)
Then the virtual power principle (2.62) becomes :

δv i M i j v j + f int i -f ext i = 0 (2.72)
Because the virtual velocity can be chosen arbitrarily, one can obtain the momentum equation for the node i :

M i j v j + f int i + f ext i = 0 (2.73)
and then the system of equations for all the nodes writes :

Ma s + f int -f ext = 0 (2.74) with M =              M 11 1 d M 12 1 d • • • M 1i 1 d • • • M 1n 1 d M 21 1 d M 22 1 d • • • M 2i 1 d • • • M 2n 1 d . . . . . . . . . . . . . . . M i1 1 d M i2 1 d • • • M ii 1 d • • • M in 1 d . . . . . . . . . . . . . . . M n1 1 d M n2 1 d • • • M ni 1 d • • • M nn 1 d              f int =              f int 1 f int 2 . . . f int i . . . f int n              f ext =             f ext 1 f ext 2 . . . f ext i . . . f ext n             (2.75)
where 1 d denotes the unit matrix of the dimension of the problem, and a s the nodal acceleration matrix of dimension [nd × 1], with d being the number of dimensions of the investigated problems and n the total number of nodes :

a s = a 1 a 2 • • • a i • • • a n (2.76)
The same decomposition is made for nodal velocity v s and displacement u s .

Discretization of constitutive relation

The semi-discrete equations (2.74) is applicable to linear or non-linear cases, since the methods for calculating the internal and external forces are not specified. For example, when calculating f int , the material property may be linear or non-linear, and the deformation of the solid may be small or large. The constitutive model of the solid has been presented in section 2.1.3.2. The constitutive equation which gives the relation between the stress and the strain is elasto-plastic in this case, but only the linear elastic will be discretized in this part, and the infinitesimal strains hypothesis is considered. This assures a linear relation between σ and ε, and between ε and u. The strain from equation (2.15) can be linearised for infinitesimal displacements and rotations, and becomes the symmetrized gradient of displacement :

ε ≈ ε 1 1 2 ∇ x u + ∇ x u = ∇ sym x u (2.77)
Then, the Voigt notation will be used to simplify the Hooke's law (2.16) σ i j = C i jkl ε kl , by writing σ and ε as vectors and C as a matrix :

{σ} = [C]{ε} (
{δD} = B i δv i {ε} = B i u i (2.82)
where u i denoting the nodal displacement for the node i. With (2.78) and (2.82), the virtual internal power can be calculated by :

δP int = Ω s (B i δv i ) [C]B j u j dΩ = δv i f int i (2.83)
Hence, the internal nodal force f int i for the node i is :

f int i = Ω s B i [C]B j dΩ u j = K i j u j (2.84)
with which the momentum equation for the node i (2.74) becomes :

M i j ü j + K i j u j = f ext i (2.85)
Then the semi-discrete equations for all the nodes in matrix form are :

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

K i j = Ω s B i [C]B j dΩ (2.87)
Note that in linear elastic cases with infinitesimal strains, the stiffness matrix K can be considered to be constant. The mass matrix M is always the same, even for the non-linear cases with large strains, because the mass matrix is calculated by using equation (2.13) into equation (2.71) :

M i j = Ω s N i (X)ρ s (x,t)N j (X) dx = Ω 0 s N i (X)ρ 0 s (X)N j (X) dX (2.88)

SPH for fluid sub-domain 2.2.2.1 Presentation of the SPH method

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

The SPH method was first proposed by Gingold and Monaghan [GIN 77] and independently by Lucy [LUC 77]. The basic idea of this method is to represent the value of a function of a certain particle i by using the values of other particles which are located sufficiently near to this particle i. The concept of the SPH method starts with the integral representation of a function f (x) :

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

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

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

f (x) < f (x) >= Ω f (x )W (x -x , h) dx (2.92)
where h is the smoothing length which defines the influence area D of the smoothing kernel function W . In the SPH convention < f (x) > denotes the kernel approximation of the function f (x). As < f (x) > is only an approximation of f (x), the choice of the kernel function is very important. Some basic conditions for the kernel function W are presented in the following :

Normalization condition D W (x -x , h) dx = 1 (2.93) Symmetry condition W (x -x , h) = W (x + x , h) ∇ x W (x -x , h) = -∇ x W (x + x , h) (2.94) Delta Dirac function property lim h→0 W (x -x , h) = δ(x -x ) (2.95) Compact condition W (x -x , h) = 0 when |x -x | > κh (2.96)
where κ is a constant defining the non-zero area of the kernel function, which is called the kernel ratio.

An example of the kernel function is shown in Figure 2.3, in which i is the investigated fluid particle, and j is one of the neighbour particles. The coordinates of the fluid particle i and j are x and x respectively, and r = x-x . Finally, κh is the radius of the influence area D i of the kernel function, which is bounded by ∂D i . The kernel function chosen in this work will be Wendland C4 kernel function, for which details can be found in [SCH 06]. In the spatial integral over the influence area, the support domain D i is approximated by the sum carried out over all the neighbour particles, then (2.92) becomes :

W(r,h) i j ∂D i D i r κ h
f i j∈D i f j W i j ω j (2.97)
where, for the sake of brevity,

f i = f (x i ), W i j = W (x i -x j , h)
, and ω j is called the weight of particle j, and denotes its volume in 3-dimensional cases (area in 2-D). Equation (2.97) is called the particle approximation.

After some developments [LI 13], the expression of the gradient of a function f over a particle i can be found by [ MON 88], which assures the zero order consistency :

[∇ x f (x)] i j∈D i ( f j -f i )∇ i W i j ω j (2.98)
where ∇ i is the gradient operator with respect to x i . The interesting point of the SPH method is that the derivation is carried out on the known kernel function, as seen in equation

             ∂ ∂t χ Ω ρ f dΩ + Ω ∇ x • ρ f (v f -v) dΩ = 0 ∂ ∂t χ Ω ρ f v f dΩ + Ω ∇ x • ρ f v f ⊗ (v f -v) + p f 1 dΩ = Ω ρ f b dΩ (2.
                                                               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 j∈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 -ω i k∈∂D i W ik 2 ρ E,ik v E,ik ⊗ (v E,ik -v 0,ik ) + p E,ik 1 • 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 1 • ∇ i W i j ω j (2.100)
where v i = v(x i ,t) denotes the Lagrangian velocity of particle i and ρ i = ρ(x i ,t) the density of fluid associated to particle i. v 0,i , v 0, j and v 0,k denote the arbitrary particle velocity vectors of particles i and its neighbour j, as well as discrete wall element k. n k is the unit normal vector pointing from the fluid to the solid, and s k stands for the surface area of the wall element k. The symmetrized velocity v 0,i j is defined as v 0,i j = 1 2 (v 0,i + v 0, j ). Symmetrized density ρ E,i j , velocity v E,i j and pressure p E,i j are the intermediate fluid status obtained by solving a Riemann problem. Concerning the density ρ E,ik velocity v E,ik and pressure p E,ik , they are the intermediate fluid status obtained by solving a partial Riemann problem at the solid wall boundary. g i stands for the acceleration due to the gravity field acting on particle i.

Time integration

Governing equations of solid and fluid sub-domains have been discretized in space in previous sections. They must be integrated in time for transient computation. The physical phenomenon involved is an impact, then the problem deals with fast dynamics, thus explicit schemes will be considered for both sub-domains.

Time integration scheme for solid sub-domain

For the solid sub-domain, the Newmark time-integrator is chosen [NEW 59]. Given the linear dynamic equilibrium equation (2.74) at the instant t = t n+1 :

Ma s (t n+1 ) + Ku s (t n+1 ) = f ext (t n+1 ) (2.101)
or, written in compact form :

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

     u n+1 s = u n s + ∆tv n s + ∆t 2 2 (1 -2β)a n s + 2βa n+1 s v n+1 s = v n s + ∆t (1 -2γ)a n s + γa n+1 s (2.103)
where u n s denotes the displacement field at the instant t = t n , v n s the velocity field and a n s the acceleration field, which are already known. ∆t = t n+1t n represents the time step. β and γ are the coefficients of the Newmark scheme. Different combinations of these coefficients give the different properties of the numerical scheme [BAT 76] :

• Unconditionally unstable if γ < 1 2 (2.104) • Unconditionally stable if 1 2 ≤ γ ≤ 2β (2.105) • Conditionally stable if γ ≥ 1 2 , β < γ 2 , ∆t ≤ Ω crit ω max (2.106)
where Ω crit for undamped system is :

Ω crit = γ 2 -β -1/2
(2.107)

and ω max corresponds to the highest natural frequency of :

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

Time integration scheme for fluid sub-domain

The 2 nd order explicit Runge Kutta scheme (RK2) is used as the time integrator for the fluid sub-domain. After the spatial discretization, the equations in (2.100) are all ordinary temporal differential equations, which can also be written as :

dQ dt = H (Q) (2.109)
where Q = Q(t) depending only on time t. In the system (2.100) for the fluid particle i, Q may be the position x i , the weight ω i , the mass ρ i ω i , or the momentum

ρ i ω i v i . The operator H (Q) is a function of Q, which represents the spatial discretization part.
The RK2 scheme is often called the two-stage Runge-Kutta scheme, because there are two stages for going from t n to t n+1 . Indeed, there is an intermediate state at t = t n+1/2 for which the variables are calculated.

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

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

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

Q n+1 = Q n + ∆tH (Q n+1/2 ) (2.111)
where the function H (Q n+1/2 ) is calculated with the variable Q n+1/2 obtained in the first stage.

Because one should calculate first the state of variables at the instant t n+1/2 , this method is also called the mid-point version of the 2 nd order Runge-Kutta scheme. However, the RK2 scheme does not have a unique form. For instance, the Heun's method is also classified as a RK2 sheme. In this work, the mid-point version of RK2 will be used for time integration of the fluid sub-domain. As an explicit time integrator is used, the Courant-Friedrichs-Lewy (CFL) condition needs to be satisfied in order to have a stable numerical simulation. In this work, the applied CFL condition is the one given in [LED 10] :

∆t ≤ K CFL min i h i c i + 2 v i 2 (2.112)
where the coefficient K CFL ∈ [0, 1] is a function of the used time integrator, which is usually determined through the practical numerical experience.

FSI coupling

The coupling algorithm used for the present simulations has been developped by Li and Nuñez-Ramirez [LI 13, NU Ñ 17b, NU Ñ 17a], and is based on the work of Combescure, Gravouil and Mahjoubi [MAH 09, MAH 09]. The numerical stability of the global problem depends on the continuity conditions prescribed at the interface [COM 02]. A zerointerface energy condition is considered because, numerical coupling is willing to generate an artificial energy loss or gain. The objective is to minimize the increment of interface energy ∆W I over the time t ∈ [t n ,t n+1 ] defined as :

∆W I = t n+1 t n Γ I n s • (-p s 1) • v s + n f • (-p f 1) • v f dΓ dt (2.113)
After mathematical transformations, the kinematic and pressure conditions presented in (2.49) become :

• Pressure condition p f (x k ) = p s (x k ) = p k (2.114)
• Kinematic condition -In geometrically linear cases, the normal n k does not change with time :

n k • v s (x n+1 k ) + v s (x n k ) 2 -n k • v f (x n+1 k ) + v f (x n k ) 2 = 0 (2.115)
2. Numerical methods for droplet impact -In geometrically non-linear cases, due to the kinematic conditions at mid-step :

1 2 n n+1 k • v s (x n+1 k )s n+1 k + n n k • v s (x n k )s n k - 1 2 n n+1 k • v f (x n+1 k )s n+1 k + n n k • v f (x n k )s n k = 0 (2.116)
where x n k denotes the position of the interface element k at instant t = t n , s n k its surface area, and

n n k = n f (x n k ) = -n s (x n k
) its normal vector. Note that x k and n k are noted without n superscript for geometrically linear cases, because they are not time-dependant in these cases, unlike non-linear ones, where they may change over the time t ∈ [t n ,t n+1 ].

The rest will not be detailled here, and more informations can be found in PhD works of [LI 13, NU Ñ 17a].

Numerical simulation of droplet impact phenomenon

The aim of this section is to present the coupled system features, and determination of the most optimized model to lighten the simulation. The fluid computation being heavy, a simplification procedure is presented, then the validation procedure is performed to corroborate the relevance of this reduction. Finally, the coupled simulation is performed and analysed in terms of load and stress history.

Model features

The solid sub-domain is computed by the Finite Element Method (FEM) with the explicit dynamics code EuroPlexus R [CEA 02], which is developed jointly by the french Commissariat à l' Énergie Atomique et aux Énergies Alternatives (CEA) and the European Commission / Joint Research Center (EC/JRC). This code is suitable for highly non-linear explicit dynamics with erosion. Concerning the fluid sub-domain, ASPHODEL code is used. This in-house code developed by Andritz Hydro uses the Smoothed Particle Hydrodynamics method (SPH) and is efficient to treat free surfaces. The Fluid-Structure Interaction (FSI) is performed by the two-way coupling code developed by Nuñez-Ramirez [NU Ñ 17a], which is energy-conservative at the interface for same time-steps.

The solid body consists of a rectangular/parallelepipedic shaped steel domain, with height H and width 2H. To avoid waves reflection, H > c L T c , where T c is the elapsed time after impact, and c L the compression waves velocity defined in equation (1.1). The material considered is homogeneous, isotropic and perfectly bilinear elastoplastic with Young's modulus E 0 , second tangent modulus E t , Poisson's ratio ν, initial yield strength σ 0 Y and density ρ 0 s , and only isotropic hardening is considered (see section 2.1.3.2 for details). The water droplet is a full disc/sphere with a radius R = 0.5 mm and moves perpendicularly towards the solid with an impact velocity V = 100 m.s -1 . The fluid core has a density ρ 0 f and an ambient sound velocity c 0 f , and is an inviscid fluid, i.e without viscous effects. Surface tension is not taken into account either (see section 2.1.4). All of the material data are summed up in Table 2.1 and the model is represented in Figure 2.4.

ρ 0 s 7700 kg.m -3 ρ 0 f 1000 kg.m -3 E 0 200 GPa c 0 f 1500 m.s -1 E t 20 GPa γ f 7 σ 0 Y 560 MPa ν 0.288

Convergence of fluid computation

In this section, the numerical convergence of the numerical model is checked, to get the most optimized discretization to reduce the model size, thus the CPU time. The maximum SPH particle size is investigated, as the possibility to consider a two-dimensional problem instead of a three-dimensional one, and an eventual vertical truncation of the fluid domain.

It is important to note that axisymmetric formulations do not exist yet for the SPH. Thus the reduction from a three-dimensional model to a a two-dimensional one is checked, not axisymmetric. In this case, the solid domain is modelled by plane strain.

Impulse as reference quantity

The numerical convergence investigation has to be build based on a reference physical quantity. This test is produced considering the impact of the water droplet on a rigid wall, representing the solid. The contact pressure p(x,t) could be compared between the different models over time and space and could work as the reference quantity. However, in fast transient dynamics, i.e for impact physics, the force F(x,t) (or pressure) is less representative than the impulse, which is the time-integral of the force :

J (x, T F ) = t 1 t 0 F(x,t) dt T F = t 1 -t 0 (2.117)
where t 0 and t 1 stand for the initial and final time respectively, T F is the time duration of the considered impulse J (x, T F ). For a constant force over time F(x), the impulse becomes simply : J (x, T F ) = F(x)T F . Moreover, the contact pressure shows high perturbations when measured, which makes the comparison between models more complex.

The impulse being the time-integral, "absorbs" these numerical flaws over time, which provides an more elegant and proper comparison.

Let us consider a one-dimensional dynamic system consisting of a mass M 1 , a damping C 1 and a stiffness K 1 , and under a time-dependant load F(t). The coordinate of the mass is noted z(t), and its velocity and acceleration ż(t) and z(t) respectively. For fast dynamics, the damping can be neglected C 1 = 0, and unit mass M 1 = 1 and stiffness K 1 = 1 will be considered for simplification.

The mass is at rest at initial time, and is loaded by a constant force F 0 at t = 0 during a time T F , which then vanishes at t ≥ T F . This system is illustrated in Figure 2.5.

The problem is solved in two distinct time phases :

z(t) = z -(t), t < T F z(t) = z + (t), t ≥ T F (2.118)
For t < T F , the equations of motion are, with the initial conditions :

       z-(t) + z -(t) = F 0 z -(0) = 0 ż-(0) = 0 (2.119)
The solution of (2.119) is :

F(t) z(t) M 1 = 1 K 1 = 1 C 1 = 0 F(t) t F 0 0 T F Figure 2
.5: Left : Schematic representation of the dynamic system. Right : Time evolution of load force.

z -(t) = F 0 [1 -cos (t)] ż-(t) = F 0 sin (t) (2.120)
When T F → 0, a Taylor series expansion with second order truncation error leads (2.120) to :

     z -(t) = 1 2 F 0 t 2 ż-(t) = F 0 t (2.121)
Now for t ≥ T F , i.e. after the impact, the equations of motion are, with respect to the continuity conditions:

       z+ (t) + z + (t) = 0 z + (T F ) = z -(T F ) ż+ (T F ) = ż-(T F ) (2.122)
Using the same method as before, for t ≥ T F , the solution of (2.122) becomes :

         z + (t) = F 0 T F sin(t) - T F 2 cos(t) ż+ (t) = F 0 T F cos(t) + T F 2 sin(t) (2.123)
The solution (2.123) is the dynamic response of the oscillator to an impact. This equation shows that the amplitude of the response is proportional to F 0 T F , which is simply the impulse of the constant force F 0 during the time T F . For a deformable body, the displacement leads to strains, leading to stresses, and the velocity to strain rate. Finally, the impulse is the quantity driving the response of the structure for impact dynamics, and the impact duration T F is a predominant factor : a higher duration generates a larger amplitude. This last assertion is true for a sufficiently small time, i.e T F 2π/ω 1 , where

ω 1 = È K 1 /M 1 = 1
is the natural pulsation of the oscillator, otherwise the problem does not deal with impact dynamics.

This study allows choosing the contact impulse as the reference quantity for the numerical convergence study of the fluid computation.

Impulse calculation

The impulse will be calculated with probes i distributed on the wall. Each of them provides the pressure over time p(x i ,t), which is then integrated over time T F to get the impulse per unit area J A (x i , T F ) :

J A (x i , T F ) = T F 0 p(x i ,t) dt (2.124)
which, for times t j with t 0 = 0 and t N = T F gives numerically by trapezoidal rule :

J A,i (T F ) 1 2 N j=0 p j i + p j+1 i t j+1 -t j (2.125)
where, for the sake of brevity, p j i = p(x i ,t j ) and J A,i (T F ) = J A (x i , T F ).

Reducing the problem dimension

The physical problem involves fast dynamics with explicit schemes. Thus, small time steps are used to ensure the numerical stability. Moreover, the jetting phenomenon requires an accurate space discretization, which reduces even more the time step. A three dimensional computation needs a large number of particles, where a two dimensional one should be less heavy. For example, considering a spheric shaped domain, the ratio between the number of particles in 3-D (sphere) and 2-D (disc) is about n R , where n R = R/∆x is the number of particles of size ∆x standing in the sphere radius R . If there are about 1000 particles inside the sphere radius, the 3-D model contains 1000 times more particles than the 2-D one. Reducing the size of the model by considering a 2-D model instead of a 3-D one should be interesting thereby.

In two dimensional cases, probes are placed on the wall on (O, x) axis. They are close to each other on the center of impact area near x = 0, and around the jetting zone x = R jet , in order to accurately capture the pressure peak signal. This is explained in Figure 2.6.

For three dimensional cases, the procedure is quite similar, but the probes are placed on three different axis of the (x, O, y) plane. Indeed, we have no certainty that the fluid behaviour would be identical in every direction and jetting could not be an axisymmetric phenomenon. The three directions chosen are : on the (O, y) axis, and the bisectors of ( x, O, y) and ( y, O,x), as shown in Figure 2.7.

Each of the three lines contains the same number of probes, placed at the same radial location from (O, z). For this reason, the mean value between the three lines is calculated for each probe from the same radial location. The results are presented in Figure 2.8. The mean value of impulse along the three lines is compared to the one calculated for 2-D model. The result is displayed in Figure 2.9.

These two results show a slight scattering of the results in 3-D model between the different directions. However, the mean value is of the same order of magnitude as the result of 2-D model. This approximation is sufficient for the present work, because the main phenomena like jetting are present in 2-D model too. The following simulations shall be performed considering a two dimensional model.

Symmetry

The problem seems to be symmetric at first view, because the geometry, the material, the loading are symmetric, and so the behaviour should be so. Considering a symmetric model would allow us to reduce the number of particles by two. In SPH method, the symmetry condition needs virtual particles, to ensure a good numerical integration for the kernel function. These particles are placed by mirror reflection, on the other side of the symmetric axis/plane. The SPH code used here manages the symmetry, but the equivalence between a full and symmetric models must be checked for the current physical problem.

For the symmetric model, the probes are placed on the positive part of the (O, x) axis. For the full one, they are placed by symmetry on the negative part too, and the value of pressure is averaged. The results are presented in Figure 2.10.

The two results are similar, thus the use of symmetric model is relevant.

Truncation

As the problem deals with fast transient dynamics, the behaviour of the physical domains is driven by the propagating waves. Inside the droplet, the compression wave starts from the contact zone at z = 0 and moves on toward opposite side of the droplet at z = 2R (cf. section 1.1.1.1). The distance travelled on the vertical axis Z by the wave front depends on the time elapsed after impact T and is Z = c f T , for a constant sound velocity (section 1.1.2.1 and equation (1.5)). With an initial sound velocity of c 0 f = 1500 m.s -1 , an impact velocity of V = 100 m.s -1 and a droplet radius R = 0.5 mm, the ratio Z/R is equal to 0.5 and 1 for a time T of 150 ns and 300 ns respectively. In other words, the wave front reaches Z = R/2 at 150 ns and Z = R at 300 ns.

If the calculation time T c is lower than the time T for which the front waves reach the heigh Z = c f T , the higher part of the droplet has no influence in fast dynamics, i.e the part z > c f T c . From this statement, it would be relevant to ignore the part of the droplet, for which z > Z, this shall reduce the number of particles again.

The truncated models at z = R are z = R/2 are compared to the full model for a calculation time of T c = 100 ns. As the wave front reaches z = R/2 at t = 150 ns, the result should not be affected. The three models are simplified by symmetry, as this this feature is relevant. The three models are presented in Figure 2.11.

The results are presented in Figure 2.12, and the impulse distribution is similar for the three models. The assumption of non-influence of the higher part of droplet is verified and the fluid sub-domain can be truncated at z > c f T c . One should note that the truncation height depends on the computation time T c , but also on the impact velocity V , since the sound velocity depends on it :

c f = c 0 f + k f V .

Particle size

Here we investigate the particle size, to find the maximum value for which the results remains unchanged. The previous computations have been performed with a size of ∆x = 6 • 10 -7 m. Other results not presented here show that a particle size in the order of ∆x = 10 -5 m gave wrong results, and in the order of ∆x = 10 -6 m gave strange results. Remark : This comparison shows the impulse at t = 64 ns after impact, because the model with ∆x 2 showed instabilities afterwards.

The impulse is similar for the two particle sizes, showing that a size of ∆x = 6 • 10 -7 m is sufficient to ensure the required numerical accuracy for the current physical problem.

Conclusion on fluid model reduction

Finally, these results show that the numerical model of the droplet can be set in two dimensions, with a symmetry on (O, z), a top truncation depending on the simulation time, and a particle size of 6 • 10 -7 m. The comparison of CPU time for the different models is given in Table 2.2 for the 2-D models only. Indeed, the 3-D model was computed on a different architecture due to its weight, and the comparison of its corresponding CPU time with others ones may not be relevant.

One-way coupling

As observed in the previous section 2.3.2, the transient analysis of droplet impact can be performed by a two dimensional model, with symmetry, a truncature depending on the physical time, and a particle size of ∆x = 6 • 10 -7 m. The details of the model are given in section 2.3.1. The structural sub-domain is described by plane strains and uses the symmetry property at x = 0, like the fluid sub-domain. The mesh is composed of quadrangles with one Gauss point, because this shape is more suited than triangles for impacts, and a single integration point has proved sufficient for the current application.

Procedure

A first computation without coupling is performed in two steps :

• Step 1 : Fluid computation. The droplet impact is performed on a rigid wall in order to record the pressure distribution history p(x,t). Numerically, this results in the pressure p j i = p(x i ,t j ) for each probe i at position x i and time t j . This time, each probe i is equidistant from its neighbours i -1 and i + 1.

• Step 2 : Structural computation. The pressure p j i is then used as load on the interface, i.e on the edge of the contact element. The normal stress of the edge of each element is equals to the pressure acting on the fluid by the rigid wall σ n (x i ,t j ) = -p(x i ,t j ). The solid mesh elements should measure the same size as the fluid particles.

Results

The pressure distribution over time is given in Figure 2.14. The pressure peak is found at x/R ≈ 0.18 and corresponds to the results found in literature (cf. section 1.1.2.2). The peak has a value of p max ≈ 1.7 GPa, which equals ten times the water-hammer pressure calculated for this impact velocity (section 1.1.2.1 and Figure 1.5). The jetting time seems to be T jet ≈ 70 ns after impact, which corroborates the fast dynamics aspect of the physical problem. The Figure 2.15 shows the pressure field inside the droplet and the jetting phenomenon. This frame corresponds to t = 150 ns after impact, so after the jetting formation, and this phenomenon is clearly apparent, and the high deformation of the fluid region shows the strength of the current SPH code.

Figure 2.16 shows the propagating compression wave inside the solid at t = 600 ns after impact. For convenience, the elements have a size of 5 • 10 -5 m. Indeed, ten probes have been selected in the radius x ∈ [0, R], thus the corresponding elements have a size of R/10.

Limitations of this method

This procedure is complex to perform because the pressure has to be recorded and transmitted in the FEM code with high time and space accuracy. The example shown here uses only ten probes, whereas there are about 10 4 particles in the droplet radius. Moreover the location of probes stops at x = R, while the interface is larger. Then, the pressure is recorded on a rigid wall, which do not fit with the real deformable solid. In case of large deformations, these pressure values would be wrong. Next, the erosion procedure requires a second computation after the first erosion mechanism (the procedure is described in Chapter 3). The pressure field should be very different consequently and the geometry of the construction new wall should be laborious.

For all these reasons, the complete analysis requires a two-way coupling, what is performed in the next section.

Two-way coupling

Using equation (1.1) with material data from Table 2.1 allows obtaining numerical values for velocities c L ≈ 5820 m.s -1 and c T ≈ 3175 m.s -1 , which are both higher than the ambient water sound velocity considered c 0 f = 1500 m.s -1 . Propagation of compression waves is given by observing the evolution and distribution of pressure inside the droplet p f . A negative pressure corresponds to tension, and a positive one to compression. Concerning the solid, hydrostatic stress σ H defined in equation (2.21) can be used . Unlike p f the sign convention of σ H is the following : σ H > 0 in regions subjected to traction and σ H < 0 to compression. Indeed, equation (2.29) gives a proportional relation between hydrostatic stress σ H and dilatation Θ s , with K s > 0. Thus, hydrostatic stress describes the traction-compression state. For solids, shear waves act jointly to compression waves (cf. 1.1.1.2). They can be observed with the Von Mises stress σ VM defined in equation (2.24). Indeed, shear is related to distortion, and σ VM depends on the stress deviator s which defines distortion energy U d in equation (2.23).

As Von Mises stress stands for the shear intensity and hydrostatic stress gives information about the straction-compression state, the signed Von Mises stress defined in equation (2.25) gives both informations about shear and traction-compression and is useful to locate zones carrying fatigue load and potential crack initiation and propagation. Indeed high shear can lead to cracks initiation and a traction state tends to open them (and a compression state to close them).

Figure 2.17 shows both hydrostatic and Von Mises stresses acting in the structure for three different times after impact. By observing those quantities, P-waves and S-waves front appear to show the velocity difference. However, the front gap between P and S-waves is not clear, due to the nature of load. Indeed, droplet impingement does not act like a proper impact, as seen in Figure 2.14, the pressure depends on time and especially on space too. Then, its maximum value is reached when microjets form, thus on contact boundary. Finally, contact area growths faster than compressions waves inside fluid during the first moments of the impact, which generates other fast moving loading sources. This complex load induces mixed moving stress, such as traction/compression and shear at the same time, and the location of loading sources moves fast. Nevertheless, the travelling of a traction-compression wave is observed, which shows a general overview of the zones carrying high stresses in term of fatigue. Analysis of hydrostatic and Von Mises stresses allows globally estimating the type of load inside the solid volume. For a solid region near the surface, the different times of Figure 2.17 show a change of sign for hydrostatic stress. Indeed, hydrostatic stress in Figure 2.18 shows that this region is subject to a traction-compression cycle during the droplet impingement. Signed Von Mises stress defined in equation (2.25) is given in Figure 2.20 at the time of jetting. Its magnitude gives the shear intensity and its sign shows the traction or compression state. Zones where high discontinuities stand are subject to combined effects of both traction/compression and shear, which are important quantities for crack initiation and propagation.

The study of impulse on section 2.3.2 shows the most deformed zones on the interface layer, since the displacement is proportional to the impulse (see equation (2.123)). The Figure 2.12 shows a maximum at the center of the contact zone (x = 0). Despite the pressure peak, which locate at x/R ≈ 0.18, the maximum shock intensity is found at a different place. Signal disturbance is observable around x/R = 0.18, because the microjets lead to local instabilities.

Results show absence of plasticity due to low stress intensity compared to the yield stress (Figures 2.17, 2.18 and 2.20). This assertion is well checked by observing the cumulated plastic strain, which is zero. Low stresses and absence of plasticity lead to high cycle fatigue domain (plasticity and hardening leading to oligocyclic fatigue domain). This information is important for the choice of fatigue criterion.

The principal stresses (presented in 2.1.3.4), which denote the maximum stress state are shown in Figure 2.21. Their evolutions corroborates the previous results, because the maximum stress reaches only 200 MPa.

The triaxiality of stresses ζ (presented in 2.1.3.3) is useful too for stress analysis. Its time evolution is presented in Figure 2.19, and the values show that both load types (traction/compression and shear) are present at the same level, because ζ does not stabilize at 0 (pure shear) or ±∞ (pure hydrostatic pressure). Thus the stress state is a combination of all types of loads.

The next features that can be observed are the principal stresses ratios r i, j defined in equation (2.40). They are plotted over time in Figure 2.22, and it appears that their evolution is not constant, particularly for r 1,2 . This means that the simulation involves non-proportional or non-radial load. This information is also predominant for the choice of fatigue criterion in Chapter 3. Proportional and non-proportional load are described in section 3.1.5. The last important quantity is the angle between the principal basis B p = (e I , e II ) and the reference one B r = (e 1 , e 2 ), calculated easily in plains strains/stresses as :

α p = cos -1 (e I • e 1 ) (2.126)
Its time evolution is plotted in Figure 2.23 and allows quantifying the rotation of principal axis. The angle is close to zero until t = 250 ns, but stresses are present (as seen in Figures 2.18 and 2.21), which can correspond to the propagation of P-waves, because traction-compression does not change the principal directions. After that, the angle starts to increase and oscillates between -90 • and 90 • , which denotes total rotation principal axis. The first high peak at t ≈ 280 ns can represent the S-waves propagation, because shear rotates the principal axis. This last analysis is important too for the definition of fatigue criterion.

Conclusion on numerical methods for droplet impact

In this chapter, the numerical methods build for droplet impact simulation have been presented. The physics is based on continuum mechanics, thus the equations of solid and fluid mechanics have been presented, such as conservation laws and constitutive/state equations. Coupling conditions needed to be added to simulate correctly the interaction between the two sub-domains. Then, the numerical disretization methods have been described. The solid is discretized by the Finite Element Method and the fluid by the Smoothed Particle Hydrodynamics. To ensure a good compromise between accuracy and CPU time saving, a numerical convergence study is performed for the fluid computation. The quantity used to evaluate the validity of results was the pressure impulse, because it drives the response of the mechanical systems in fast transient dynamics. Results have shown that considering a two dimensional problem was relevant. One only considers a symmetric model, and a truncation of the top on the droplet, depending on the simulation time and the impact velocity. The particle size has been investigated too, to find the best equilibrium.

Next, a one-way computation has been build in two steps. A first simulation with the droplet impacting a rigid wall was performed, and the time distribution of contact pressure on wall was recorded. The second simulation involves the structure, and the previous pressure is entered as a load on the interface edges of the elements. This procedure allows understanding the general behaviour the whole system, but is complex to implement, and not adaptable to modifications of interface shape. On this basis, considering a two-way coupling should be preferred. This last simulation is a stronger model for analysing the structure behaviour. The stresses involved are not high enough to cause a sudden fracture of the material. Nevertheless, the regions near the contact are loaded by alternating traction and compression, which is a predominant factor for a fatigue cracking.

This simulation allowed corroborating the fatigue cracking as damage mechanism for the current droplet impact phenomenon. The next chapter will present the implementation of a fatigue analysis in the numerical model. Chapter 3

Fatigue developments

This chapter presents the necessary elements used for the fatigue analysis, which is the damage mechanism of the current study. First, the fatigue phenomenon is presented, then a fatigue criterion is selected to define the eroded/non-eroded state of material, and finally we explain the implementation procedure of the fatigue analysis into the simulations.

Generalities

Fatigue failure consists in a continuous succession of several complex phenomena : hardening, local damage of the material and leads to the birth of first microcracks (incubation phase), then the slow or fast advance of these microcracks resulting in fracture of material (propagation phase). These steps correspond to the ones defined in Chapter 1, section 1.2.2.1, for erosion by water droplet impingement, since the erosion mechanism is fatigue cracking.

The local nature of fatigue damage requires the study of load paths, and the domains with high stresses state must be highlighted. Most metallic materials are polycrystals. Thus, the diversity of grains orientations, dislocations, precipitates, second phases, and other micro-irregularities generate local anisotropy and heterogeneity inside the material. The fatigue damage starts from grains or precipitates which are the most unfavourably oriented to the maximum stresses or strains directions. Then, these cracks propagate on particular crystallographic planes. This is the reason why metallic fatigue is strongly linked to the material microstructure.

More details about the fatigue cracking may be found in [BRA 84, LU 02].

Definition of terms used in fatigue

Generally, fatigue experimental tests are performed on a test coupon, which is loaded in traction-compression, shear, torsion or bending, or combinations of them. Normal stresses are generally denoted by σ, and tangential ones by τ.

The sample is then loaded, with several stress cycles, i.e the stress, with a variable amplitude repeated periodically. The most used is the sinusoidal cycle. The description 3. Fatigue developments below uses the normal stress, but works with tangential stress too.

For a given load, the stress can be expressed as a combination of two terms : the mean stress σ m (constant) and the alternate stress σ a (time-dependant) :

σ(t) = σ a (t) + σ m (3.1)
The mean stress is for a sinusoidal cycle the mean of maximum and minimum values, respectively σ max and σ min :

σ m = σ max + σ min 2 (3.2)
and the alternate stress is equal to :

σ a (t) = σ am sin(ω σ t) (3.3)
where σ am is the stress amplitude, and ω σ denotes the signal pulse of the stress load. The stress amplitude is then defined as :

σ am = σ max -σ min 2 (3.4)
The stress range ∆σ = 2σ am is sometimes used instead. All these quantities are shown in Each of quantities σ max , σ min , or σ m can be positive, negative or zero 1 . To distinguish all these cases, the load ratio R σ is defined as the ratio between the minimum and The fatigue phenomenon maximum stresses :

R σ = σ min σ max (3.5)
Three extreme cases are described :

• R σ = -1 corresponds to a perfect symmetric alternate cycle, because σ m = 0.

• R σ = 0 denotes a pure traction alternate case, where σ ≥ 0.

• R σ = 1 is equivalent to a constant stress or non-alternating cycle with σ am = 0.

Tests are performed for given mean stress and stress amplitude σ m , σ am until fracture. The number of corresponding cycles is N, called number of cycles before failure. The results are reported on a graphic with σ am as a function of N, called S-N diagram, or Wöhler line, illustrated in Figure 3.2.

Those tests allow determining two quantities, which describe the resistance to fatigue of a material for a given mean stress :

• Endurance limit : the highest stress amplitude that a material can survive for an infinite number of load cycles. It is noted σ D when σ m = 0, and called pure alternate endurance limit, and σ aD when σ m = 0, and called alternate or repeated endurance limit. For steels, experimentations show that endurance limit can be set to 10 7 cycles, but for high frequency in gigacyclic domain (this term is defined in next section), this limit does not exist [BAT 99]. Non-ferrous metals, like aluminium do not show an endurance limit.

• Endurance strength : the stress amplitude that a material can survive for a given number of load cycles N. It is then noted σ D (N). Concretely, endurance strength is the stress amplitude for which the material whose damage is 0.5 for a given N. For steels, endurance strength for 10 7 cycles can be mixed with endurance limit.

These quantities are different for traction and torsion. One defines σ -1 and τ -1 as the endurance limit under symmetrical alternate traction and torsion respectively. The subsript "-1" means R σ = -1 and corresponds actually to symmetrical alternate load. As measures for τ -1 are difficult to provide, the approximative value τ -1 = σ -1 / √ 3 is considered in this thesis, because it establishes an adequate approximation for most of stainless steels.

Informations about fatigue test can be found in [RAB 00a, RAB 00b].

The different fatigue domains

The mechanical fatigue can occur at any stress level. There are generally two main modes depending on the stress level according to the yield strength of the material [LEM 94, LEM 09] :

• High number of cycles fatigue : with σ < σ 0 Y , the strains are elastic or quasielastic, and the corresponding number of cycles before failure is higher than ≈ 50 000, and can reach 10 6 , 10 7 or even 10 9 , which corresponds to Gigacyclic fatigue. Depending on the field of application, the lifetime can be considered as "infinite", i.e. the loads are lower than a given fatigue limit.

• Low number of cycles fatigue : where σ > σ 0 Y , called also oligocyclic fatigue, and corresponds to numbers of cycles lower than ≈ 50 000, for which plastic strains occur too. At high temperature, fatigue phenomenon is coupled with creep damaging and can be isolated only by performing experiments at sufficiently high frequency (> 5 -10 Hz). The results of transient simulation in Chapter 2, section 2.3.4, show that stresses are far lower than the yield strength of the material. Thus our problem deals with high-cycle fatigue.

This classification is shown in

The different fatigue loads

For 1D models generally used in experimental investigations to build S-N diagrams, the stress tensor σ contains only one component σ. So minimum stress, maximum stress, mean stress and so stress amplitude are directly calculable, like in section 3.1.3. But the current case deals with multiaxial case. Indeed, we saw in Chapter 2, section 2.3.4, that all principal stresses have the same order of magnitude, which means that the fatigue analysis must consider several stress components, to calculate an equivalent stress.

Then, the type of load has its influence on how to tackle the fatigue analysis. Loads can be splitted into two types :

• Proportional load : When the principal stresses ratios r 1,2 , r 2,3 and r 3,1 remain constant over time, and principal axis do not rotate with time, the load is called proportional or radial. A typical example of proportional load is biaxial traction in phase, or traction-torsion in phase.

• Non-proportional load : In other cases, for instance if loads are sinusoidal out of phase or random. The principal axis may also rotate over time, and the load is called non-proportionnal or non-radial. This case can make stress analysis more complicated, because principal directions must be computed at any time. This type of load can be obtained for example with a test of torsion-bending out of phase combination.

As seen in Chapter 2, section 2.3.4, the ratios between principal stresses are not constant at all over time (cf. Figure 2.22), thus the load is non-proportional. Moreover, the principal axis rotate highly with time (cf. 2.23).

Conclusion on fatigue phenomenon

As seen in previous sections, our problem deals with multiaxial, high-cycle fatigue, with a non-proportional load. These three features allow defining the strategy for the fatigue analysis and the choice of fatigue criterion producing an equivalent stress, which can be used in S-N diagram to calculate the lifetime of the material.

Multiaxial fatigue criteria

General ASME criterion

One of the first historical criteria implemented in numerical codes is the one from the American Society of Mechanical Engineers (ASME) [ MAL 08]. It is based on the maximum principal stress in absolute value (spectral radius of stress tensor), like Rankine criterion but applied to fatigue. The use of principal stresses is justified because it gives the maximum stress state inside the material at a given point.

Let us consider two physical times t m and t n , with (t m ,t n ) ∈ [0, T c ], where T c denotes the largest time of computation. The corresponding Cauchy stress tensors to t m and t n are noted σ (m) and σ (n) respectively. One defines the difference stress tensor ∆σ (m,n) for the couple (t m ,t n ) :

∀ (m, n) ≤ n T , m = n ∆σ (m,n) = σ (m) -σ (n) (3.6)
where n T denotes the number of numerical time steps considered for the computation of equivalent stress. The Rankine criterion is then applied to ∆σ (m,n) in order to get the maximum value of the stress state, which is the stress range ∆σ (m,n) between times t m and t n :

3. Fatigue developments

∆σ (m,n) = Rankine ∆σ (m,n) (3.7)
This requires to diagonalize ∆σ (m,n) and find its eigenvalues ∆σ (m,n) K

, with K = I, ..., III, which denotes the principal difference stresses. Equation (3.7) becomes then :

∆σ (m,n) = max K § |∆σ (m,n) K | ª (3.8)
The maximum range over the load cycle ∆σ is then obtained by sweeping all the combinations of couples (t m ,t n ) :

∆σ = max m,n m =n ∆σ (m,n) (3.9)
Finally, the stress amplitude σ am AS equals the half of the stress range ∆σ and is :

σ am AS = 1 2 max m,n m =n § max K § |∆σ (m,n) K | ªª (3.10)
Remark : Note that in equation (3.10), the two max{•} could be swapped. The order presented here is the most natural, but is less efficient in term of CPU time. Indeed, the sweep on K requires 3 comparisons, and the sweep on (m, n) requires n T (n T -1)/2 comparisons. With the order written here, the total number of comparisons is 3n 2 T (n T -1)/2, whereas it is only 9n T (n T -1)/2 for the second solution.

This criterion shall be kept in this study because of its ease of implementation, and will be compared to another criterion provided for high cycle fatigue, described in the next section.

Specific criteria

In this section, more specific criteria are presented. Their can be structured into three main families, according to their assumptions or physical quantities involved. The description provided here do not go into details, which can be found in [CHA 05, LU 02, WEB 99b].

Empirical criteria

These criteria are build from experimental results by a general consideration of a particular multiaxial load such as traction-torsion or bending-torsion, in phase most of the time. These criteria give relevant results for particular multiaxial load, but are restricted to these specific cases only, because the stress formalism is directly linked to the type of load involved. Thus, they are not relevant in case of a general multiaxial problem.

Critical plane criteria

These criteria consider that damage is located on a certain plane, called critical plane, which depends on the current multiaxial cycle. The critical plane is supposed to contain the crack birth, and its definition is based on the choice of parameters from stresses acting on this plane, and parameters build with certain invariants associated with the stress tensor or its deviator. The choice made by authors of these criteria are based on experimental observations of fatigue behaviour of metallic materials [WEB 99b].

The terms used by these criteria are mainly the normal and tangential components of stresses acting on a plane of normal h over the multiaxial cycle considered. They are obtained by projecting the stress states onto the plane or its normal (Figure 3.3). At time t, the stress tensor σ(t) acts on the plane of normal h by its stress vector t h (t) defined as :

t h (t) = σ(t) • h (3.11)
which is decomposed into its normal component σ h (t), called the normal stress on h and its tangential one τ h (t) worn by u h and v h , called the shear vector, where (u h , v h , h) is a direct basis linked to the plane 2 :

σ h (t) = h • t h (t) τ h (t) = t h (t) -σ h (t)h = τ uh (t)u h + τ vh (t)v h (3.12)
where τ uh (t) and τ vh (t) are the respective component of τ h (t) on u h and v h . The euclidean norm of τ h (t) is called the shear intensity τ h (t) and is :

τ h (t) = τ h (t) 2 = τ 2 uh (t) + τ 2 vh (t) (3.13)
Figure 3.3 provides an illustration of the decomposition of stress vector t h (t) onto physical plane (π) of normal h at point P at time t.

Terms related to normal stress

The different criteria consider several quantities built from the normal stress σ h (t) over the cycle. These quantities are expressed in equation (3.14) :

                                 σ max h = max t∈[0,T c ] {σ h (t)} σ min h = min t∈[0,T c ] {σ h (t)} σ m h = σ max h + σ min h 2 σ am h = σ max h -σ min h 2 σ a h (t) = σ h (t) -σ m h (3.14)
2 Vectors u h and v h have no connection with displacement u and velocity v from Chapter 2.

97

P h t h (t) ¡ h (t) h (t) ¢ u h v h £ uh (t)
£ vh (t)

Figure 3.3: Decomposition of stress vector t h (t) in basis (u h , v h , h) linked to physical plane (π) of normal h. [WEB 99b]
where, σ max h , σ min h and σ m h denote respectively the maximum, minimum and mean normal stresses, σ am h is the normal stress amplitude, and σ a h (t) stands for the alternate part of normal stress at time t (or simply alternate normal stress). We can note that σ am h = max t∈[0,T c ] {σ a h (t)}. Figure 3.4 provides an overview of these quantities depending on the time-evolution of σ h (t).

σ h t σ h m σ h min σ h max period σ h am σ h am σ h a (t)

98

Terms related to tangential stresses Terms associated to evolution of tangential stress τ h (t) have a similar definition to the ones dealing with normal stress. However, their determination is more complex because of the two-dimensional characteristic of tangential stress (see equations (3.12) and (3.13)). In the same way as for normal stress with equation (3.14), the alternate part of tangential stress τ a h (t) is defined as :

τ a h (t) = τ h (t) -τ m h (3.15)
where τ m h denotes the mean tangential stress. The tangential stress amplitude is the the maximum value of the norm of alternate tangential stress τ a h (t) = τ a h (t) 2 :

τ am h = max t∈[0,T c ] {τ a h (t)} (3.16)
From a geometric point of view, calculation of stress amplitude τ am h and mean stress τ m h , requires building the smallest circle C circumscribed to the load path T in the tangential stresses space [WEB 99a], which is formed by the end of the shear stress vector τ h (t) on the physical plane over the cycle. On the contrary, there are other criteria, which maximize a damage indicator for each plane, or investigate the maximum values of stress terms. This ensures finding the plane which shows the maximum damage, in the sense that is is defined by the criterion. This is the case of the following criteria : Yokobori

Energy criteria

This last criteria family is mainly based on the invariants of stress tensor, or its deviator (cf. section 2.1. This family of criteria deals with various models and situations, because their are based on multiple theories and assumptions, however their implementation is complex. Moreover, their are designed for mainly equi-damaged planes and could not be relevant for other cases, where only some planes are activated.

Choice of fatigue criterion

The choice for our physical problem is the Dang Van criterion [DAN 71, DAN 84, DAN 89], because it is largely used in the industry, designed for high-cycle fatigue, and works with proportional or non-proportional loads and so for fix or mobile principal axis. 100 This criterion is built from microscopic estimations, considers macroscopic quantities, and is based on the following experimental considerations :

• Shear of the most unfavourably oriented crystallographic planes is the basis mechanism responsible for fatigue crack initiation, hence the importance of linking the maximum local (microscopic) shear τ a h (t) produced by the load (macroscopic) shear τ h (t).

• Hydrostatic stress σ H (t) is predominant and allows corroborating experimental results.

The criterion can be expressed by a linear relationship :

τ a h + α DV σ H = β DV (3.17)
where α DV and β DV are material dependant constants :

α DV = 3 τ -1 σ -1 - 1 2
β DV = τ -1 (3.18)
To use the criterion, a validity condition α DV > 0 is set for the material, which is equivalent to :

τ -1 σ -1 > 1 2 (3.19)
what is indeed checked for the assumption τ -1 = σ -1 / √ 3 considered previously in section 3.1.3. It can be noted that α DV no longer depends on the material.

The relation (3.17) allows defining the endurance domain which is expressed as the function of stress f DV (σ) :

f DV (σ) = τ a h + α DV σ H -β DV (3.20)
When f DV (σ) ≥ 0, the material is locally damaged on plane h. This domain is represented in Figure 3.6, where G is the load path in (σ H , τ a h ) space over time. On this graphic, two typical load paths G 1 , G 2 are represented. Both have the same shear range ∆τ h = 2τ am h and the same mean hydrostatic stress σ m H , but only G 2 causes damage. Mean hydrostatic stress is computed in the same way as normal stress from equation (3.14).

The condition of non damage can be set up as :

τ a h (t) + α DV σ H (t) β DV < 1 (3.21)
We define thereby a damage indicator E h for each plane h, as the maximization of equation (3.21) over time :

E h = max t∈[0,T c ] ¨τa h (t) + α DV σ H (t) β DV « (3.22)
The maximization of E h over h represents the determination of the critical plane, and denotes the fatigue function E DV 1 :

E DV 1 = max h {E h } (3.23)

Determination of the critical plane

Inserting (3.22) into (3.23) leads to :

E DV 1 = max h ¨max t∈[0,T c ] ¨τa h (t) + α DV σ H (t) β DV «« (3.24)
We see that there are two maximizations to perform for Dang Van criterion :

• The first one consists to find the maximum value of the term τ a h + α DV σ H . The maximum values of shear and normal/hydrostatic stress should be checked together, but it is still possible to check them individually : for each hydrostatic stress σ H (t n ) at time t n , the shear amplitude τ a h (t n ) is calculated, as seen in section 3.2.2.2. Algorithms for shear amplitude determination are proposed by [DAN 89, PAP 87, ZHO 97].

• The second one deals with checking the maximum value of the previous term for all the planes of the space. This part is the heaviest, because the angle increment must be very small during the sweep [CHA 05].

For this last reason, Dang Van [DAN 89] proposes a new method avoiding the angle sweeping. This method is described in the next section.

Second version of Dang Van criterion [DAN 89]

In the Dang Van 2 criterion, the deviatoric stress s(t) is considered instead of the shear vector τ h (t). Indeed, the deviatoric stress deals with shear stress (cf. Chapter 2). The deviatoric space being six-dimensional in 3-D, and four-dimensional in plane strains, the determination of the amplitude is equivalent to calculate the radius of the smallest hypersphere circumscribed to the load path. In this case, the load path means the set of stress states, in deviatoric space, over the current cycle. The concept of distance in the deviatoric space is set with the second invariant of deviator J 2 (cf. section 2.1.3.4). Indeed, J 2 is defined as :

J 2 = 1 2 tr (s) 2 = 1 2 i, j s 2 i j = 1 2 i = j s 2 ii + 2s 2 i j (3.25)
and the radius of the hypersphere is simply √ J 2 . Then, the terms relative to deviatoric stress are defined, in the same way as previously. The mean deviatoric stress s m is defined as the center of the smallest hypersphere circumscribed to the load path in deviatoric space. Then, the altenate deviatoric stress s a (t) is found as :

s a (t) = s(t) -s m (3.26)
The alternate shear is no longer linked to any plan h and is now noted τ a (t), without subscript. It is got by applying the Tresca criterion over the alternate deviatoric stress s a (t), and is the highest shear on principal planes (planes formed by two principal axis). Indeed, it is well known that planes of maximum shear stress are inclined at 45 • with respect to the planes of principal stress, and principal axis of s(t) are the same as σ(t) at each point (cf. section 2.1.3.4). So the alternate shear is :

τ a (t) = Tresca {s a (t)} = max K {τ a K (t)} (3.27)
where τ a K (t) are the shears on principal planes, with K = I, ...III. The classical value of shears on principal planes is given as a function of principal stresses, but according to [DAN 89], at local (or microscopic) scale, the Tresca criterion applied to stresses and deviatoric stresses is the same : Tresca {σ micro } = Tresca {s micro }, hence the use of principal alternate stresses :

τ a K (t) = |s a L (t) -s a M (t)| 2 with (K, L, M) = I, ..., III K = L = M (3.28)
where s a K (t) are the principal alternate deviatoric stresses, with K = I, ...III. They are obtained by diagonalizing s a (t) (cf. section 2.1.3.4). Alternate shears on principal planes are then :

τ a I (t) = |s a II (t) -s a III (t)| 2 τ a II (t) = |s a III (t) -s a I (t)| 2 τ a III (t) = |s a I (t) -s a II (t)| 2 (3.29)
and the developped form of alternate shear τ a (t) is :

τ a (t) = 1 2 max {|s a I (t) -s a II (t)| , |s a II (t) -s a III (t)| , |s a III (t) -s a I (t)|} (3.30)
Finally, the fatigue function for Dang Van 2 criterion contains only the maximization over time of the following equation :

E DV 2 = max t∈[0,T c ] ¨τa (t) + α DV σ H (t) β DV « (3.31)

Fatigue erosion predicting tool

General procedure

The fatigue analysis consists to a post-processing, which after the transient computation selects the eroded zones of the structure. The condition of erosion applies to elements of the mesh and is decided by the fatigue criterion.

This condition is determined by the "fatigue function" and consists to a dimensionless quantity E ∈ [0, 1], the condition of erosion being E = 1, and E = 0 corresponding to a virgin element. A fatigue criterion needs a number of load cycles N lim as input, which corresponds to lower limit for non-eroded elements. Fatigue criteria give the opportunity to predict how many identical load cycles N each element can carry before failure. If N ≤ N lim , i.e E = 1, the element is eroded. This method is a predicting one, and saves a lot of time, because only one load cycle is simulated and not N lim . After removing eroded elements of the mesh, a new FSI interface is computed and another simulation is launched for another number of cycles. The main procedure is detailed in Figure 3.7.

The resistance to a given number of cycles N lim is checked with the S-N diagram (cf. 3.2), which depends on the material. For a given stress amplitude σ am , the corresponding number of cycles to failure N is found. If N ≤ N lim , fatigue cracking initiates. Let us call σ lim the corresponding stress to N lim . The crack initiation condition N ≤ N lim is equivalent to σ am ≥ σ lim , and to E = 1. S-N curve representing a bijection, let us call the function f SN which give the stress corresponding to a number of cycles, and f -1 SN it reverse function (Figure 3.8). Obtaining the erosion condition is then performed as follows :

1. The user chooses the lower limit in number of cycles N lim 3. After computation, the stress amplitude σ am is calculated with the fatigue criterion 3. Fatigue developments for each element of the mesh for the cycle.

4. The fatigue function for each element is then computed as the ratio between the stress amplitude σ am and the stress limit σ lim :

E(N lim ) = σ am σ lim (3.32) 5.
Elements where E = 1 are considered as eroded for the current load cycle according to the number of cycles chosen N lim and the fatigue criterion. If E exceeds 1, its value is usually brought back to 1 for more relevance. Indeed, from E ≥ 1, the element concerned is damaged, regardless of the value of E.

6. The number of cycles before failure for each element can be determined as N = f -1 SN (σ am ). The fatigue criteria chosen are : General ASME criterion (section 3.2.1) and Dang Van 2 criterion (section 3.3.3). For the ASME criterion, the fatigue function E AS (N lim ) is directly given as :

E AS (N lim ) = σ am AS σ lim (3.33)
The upper stress limit σ lim is set equal to the endurance limit under symmetrical alternate traction σ -1 , and as

β DV = τ -1 = σ -1 / √ 3 (cf. 3.1.
3), we can define the Dang Van equivalent stress σ eq DV 2 and its corresponding stress amplitude σ am DV 2 :

     σ eq DV 2 = max t∈[0,T c ] {τ a (t) + α DV σ H (t)} σ am DV 2 = √ 3σ eq DV 2 (3.34)
thus its corresponding fatigue function E DV 2 (N lim ) is :

E DV 2 (N lim ) = σ am DV 2
σ lim (3.35) This equivalence of definition between E AS (N lim ) and E DV 2 (N lim ) is made for fairness and makes the programming easier and more elegant.

Protocol for stress amplitude calculation

This section presents the global overview of the different steps in the determination of the stress amplitude for each criterion.

ASME general criterion

The procedure starts with the couple of Cauchy stress tensors σ (m) and σ (n) at times t m and t n and concludes with the stress amplitude σ am AS . The protocol is described in Figure 3.9 and is performed individually for all the elements of the mesh (one Gauss point).

Dang Van 2 criterion

Here, the procedure starts with the Cauchy stress tensor σ(t) ≡ σ i j (t), and finishes with getting the stress amplitude σ am DV 2 . The protocol is presented in Figure 3.10 and is performed individually for each element of the mesh (one Gauss point).

Procedure for new FSI boundary determination

As previously mentioned, the interface between solid and fluid changes between two simulations and its shape is calculated as follows :

1. Sweep all FSI interface elements, if an interface element belongs to an eroded element, it is deleted.

2. Sweep all eroded elements and search for three sides of each element which also belongs to a non-eroded element. This will be the new FSI interface element.

This procedure is illustrated in Figure 3.11.

Conclusion on fatigue developments

In this chapter, the fatigue phenomenon is presented in a global point of view, then the main notions and quantities significations are explained. Then the different classifications of fatigue as function of stresses and type of loads are distinguished. The 2-way coupled transient computation of Chapter 2 deals with high-cycle fatigue with non-proportional load. These properties allow choosing relevant fatigue criteria : the general ASME criterion, which is built from general considerations but is easy to implement, and the Dang Van fatigue function. This analysis is useful for the current erosion simulation. Indeed the procedure consists of a post-processing computation after the FSI simulation. This operation performs a calculus on every element of the solid mesh using the previous fatigue criteria, its result allowing finding if each element is eroded or not for a given number of load cycles. And this condition of erosion is determined by the fatigue criteria. Eroded elements will be converted to "ghost", i.e they are not considered in next simulations.

This erosion simulation procedure shall be a predictive method, because only one FSI simulation is performed, and we make an extrapolation for a higher "virtual" number of load cycles.

Chapter 4 Prediction of damage due to repeated droplets impacts

This chapter presents the use of the method of fatigue analysis described in previous chapter in order to predict the damage due to droplet impact. A parametric analysis for several droplet diameters and velocities is performed, then the influence of coating layer is investigated.

Features and assumptions

This chapter presents the numerical simulations performed for prediction of the damage due to water droplet impingement. First, a parametric analysis is built for several droplet diameters φ i = 2R i and velocities V j . For each case of couple (φ i ,V j ), erosion is evaluated with both criteria ASME and Dang Van 2 and influence of each parameter on erosion is determined in terms of damage depth, shape, lifetime. The tested velocities are V = {100 m.s -1 , 200 m.s -1 } and the diameters are φ = {0.5 mm, 1 mm, 2 mm}, in order to have a ratio of two between consecutive values. The mechanical properties of the bucket material X3CrNiMo13-4 do not change from Chapter 2 (cf. Table 2.1 in section 2.3.1). Its mechanical features are given in Table 4.1, its traction curve in Then, one case is tested when adding a coating layer, which is supposed to improve the resistance of the bucket and increase its lifetime (cf. Chapter 1, section 1.2.4). The coating layer is considered as an homogeneous domain without porosity. According to ANDRITZ experts in material science, this approach should be sufficient for the present fatigue analysis. Both materials are modelled as two separated domains, without special junction law (Figure 4.1). The coating layer in composed of tungsten carbide WC-CoCr and rectangular-shaped with a thickness of 300 µm. Its spatial discretization is performed with the same elements as the raw material : 1 Gauss point quadrangles. The coating is assumed to be elastoplastic, and its material properties are presented in table 4.1 (the raw material has the same mechanical properties as in Chapter 2, cf. Table 2.1). Its traction behaviour does not assume a simple tangent modulus after plasticity, so its (σ, ε) curve is given in Figure 4

σ 0 Y 560 MPa σ 0,coat Y 1.11 GPa
The coating is supposed to be perfectly attached to the substrate, and the absence of constitutive relation for interface between raw material and coating is justified by the damage mode. Indeed, the coating is not removed by debonding or delamination, but by continuous wearing extending inside the raw material. The clean delimitation between the two domains is justified by the non-penetration of the coating in the bucket, which is actually observed in Figure 1.34 from Chapter 1.

As mentioned in Chapter 3 (section 3.4.1), the application of fatigue criteria leads to the number of cycles before failure N. This number represents the number of droplets impacts and can be expressed in terms of working hours with the following assumption : there is a droplet at every water jet passage on the current bucket, impacting on the same area. For a six-jets Pelton turbine, n jet = 6, rotating at Ω r = 400 rpm, the impact frequency 4. Prediction of damage due to repeated droplets impacts equals 40 Hz, thus the lifetime in hours L h is :

L h = 1 60 N n jet Ω r = N 144 000 (4.1)
As the stress amplitudes are noted σ am AS and σ am DV 2 for ASME and Dang Van 2 criteria, their corresponding number of cycles to failure are noted N AS and N DV 2 . The corresponding lifetimes in hours are L h AS and L h DV 2 respectively, idem for fatigue functions E AS and E DV 2 .

Results for non-coated surface 4.2.1 Results of parametric analysis

In this section, the results of the different cases (φ i ,V j ) are presented and compared to each other in terms of maximum stresses, minimum number of cycles to failure, erosion depth for a given number of impacts, and the corresponding damage shape. This analysis allows determining qualitatively the influence of each parameter φ and V in erosion intensity. Then, the difference between both criteria is investigated.

Equivalent stress and number of cycles

First, the stress amplitude σ am and the corresponding number of cycles N are computed by the fatigue program. Figure 4.3 shows the stress amplitude for both criteria (σ am AS and σ am DV 2 ) for the case (φ = 1 mm,V = 100 m.s -1 ). The corresponding number of cycles by the S-N diagram is shown in Figure 4.4. The location of higher stresses (lower number of cycles) is in this case located as the external envelope of the theoretical region of microjets birth (cf. Chapter 1, section 1.1.3).

For each case (φ i ,V j ), the maximum stress amplitude σ am is collected over the mesh for each criterion after fatigue computation. The values are shown in Figure 4.5. Figures 4.5.a and 4.5.b are equivalent, but observing them separately allows determining more clearly the influence of each parameter φ and V on the maximum stress according to each criterion. The current values of stress amplitude are always lower than the yield strenght of the material σ 0 Y = 560 MPa, which corroborates the previous assumption of high-cycle fatigue. Using the Wöhler line (S-N diagram) of the current material allows calculating the corresponding number of cycles : N = f -1 SN (σ am ). These results show that in terms of maximum stress, ASME criterion is tighter than Dang Van 2 one for all cases, and the gap increases with impact velocity V , but remains constant while droplet diameter φ changes. ASME criterion produces a higher maximum stress, but this feature is not necessary linked to the number of elements eroded for a given N lim .

The most aggressive framework is (φ = 2 mm,V = 200 m.s -1 ) (using ASME criterion), which seems logical, as kinetic energy of the droplet before impact per unit length is :

E k = 1 2 m L V 2 = π 8 ρ 0 f φ 2 V 2 (4.2)
where m L denotes the mass per unit length of the droplet1 . Thus the kinetic energy of the droplet increases with its diameter and velocity. By deductive reasoning, the most favourable case is (φ = 0.5 mm,V = 100 m.s -1 ), which is indeed confirmed by results.

The number of cycles N can be converted into working hours L h with formula (4.1). For the two extreme cases (φ = 0.5 mm,V = 100 m.s -1 ) and (φ = 2 mm,V = 200 m.s -1 ), the first element breaks after L h AS = 45 h and L h DV 2 = 230 h, and L h AS = 0.07 h and L h DV 2 = 0.47 h respectively. These last values are not relevant because the velocity of V = 200 m.s -1 is not really common. Indeed, for a Pelton height of H p = 2000 m, the corresponding water jet velocity is V p = È 2gH p = 200 m.s -1 , but considering the rotation of the runner, the relative velocity of the waterjet on the bucket is the half, so V j→p = 100 m.s -1 . Moreover, elements have a size of 10 -6 m, which represents a small erosion depth. Finally, real Pelton buckets are covered with a coating layer (the comparison is established in section 4.3), which increases their lifetime. The lifetime for the same case with V = 100 m.s -1 gives L h AS = 13.4 h and L h DV 2 = 30 h. These results show the large range of fatigue modes involved in this study : the case (φ = 2 mm,V = 200 m.s -1 ) gives an order of magnitude of 10 4 cycles which crosses the low-cycle fatigue limit whereas the case (φ = 0.5 mm,V = 100 m.s -1 ) produces 10 6-7 cycles, which almost deals with gigacyclic fatigue.

Prediction of damage due to repeated droplets impacts

N [-]

Applying erosion condition

In this section, the erosion condition N ≤ N lim is applied in order to remove eroded elements and determine the erosion intensity. The chosen lower limit will be N lim = 5 • 10 7 cycles, in order to erode the structure for all cases. The corresponding fatigue function is then determined as (cf. section 3.4.1) :

E AS (N lim ) = σ am AS σ lim E DV 2 (N lim ) = σ am DV 2 σ lim (4.3)
The erosion is characterised by the number of eroded elements, thus its depth and shape. These different aspects are analysed for all cases (φ i ,V j ).

Number of eroded elements

The total number of eroded elements for N lim = 5 • 10 7 cycles is given in Figure 4.7 for the symmetric part (the half of the solid domain, i.e on x ≥ 0). Once again, velocity 4. Prediction of damage due to repeated droplets impacts has larger influence than diameter, but both make the number of eroded elements increase. This time, the Dang Van 2 criterion seems to be harsher than ASME one. We can conclude that Dang Van 2 criterion provides lower maximum stress, but high stresses might be more spread in the solid domain, unlike ASME criterion, which produces more concentrated stresses. One can note that that the case (φ = 1 mm,V = 100 m.s -1 ) generates more eroded elements than the case with V = 200 m.s -1 for Dang Van criterion, which seems not relevant. But this result could be explained with the erosion shape associated to mesh local density. Indeed, the mesh is not entirely homogeneous and for CPU time saving considerations, its size increases when receding from the contact zone. Thus a region composed of larger elements provides less broken elements when eroded.

Erosion depth

Another important parameter is the depth of erosion, which is formed by the broken elements. It differs from the total number of eroded elements, since this number does not provide any information about the shape or dimensions of the eroded region. Figure 4.8 gives this depth for different impact velocities and droplet diameters. The same trend is observed as for previous parameters : velocity has a higher influence on the depth than diameter. Once again, ASME criterion application generates more damage, but the gap between both criteria is low compared to the depth.

Erosion shape and location

The last property that is willing to change between the different cases and criteria is the envelope of the region(s) formed by eroded elements. This is what we call here erosion shape. Previous paragraphs have presented the dimensions of the damaged region from a quantitative point of view, now the analysis focuses on the location and form of the eroded region. Figures 4.9, 4.10, 4.11, 4.12, 4.13 and 4.14 show for the different cases the fatigue function (also called damage indicator here) for the limit mentioned above N lim = 5 • 10 7 cycles, and the mesh cleared of eroded elements, i.e. elements in which E ≥ 1. We observe that the shape of erosion is more localised under the contact zone for low velocity (100 m.s -1 ), whereas for high velocity (200 m.s -1 ), it spreads also on the interface far from the contact zone.

With bot criteria, the eroded zone for low velocity is clearly enveloping the jetting (microjets) area, like a crown, which shows that microjets cause high stress intensity during their formation. The case (φ = 2 mm,V = 200 m.s -1 ) provides a circular erosion shape in the (x, O, z) plan for both criteria and seems to follow the envelope of P/S waves. Indeed, the stresses calculated for this case have a sufficiently high intensity to exceed the limit σ lim on a wide zone, i.e the region covered by the P/S waves. Both criteria generate globally the same shape for low velocity, but the difference increases for high velocity. This difference is due to the different assumptions made for each of them : the ASME criterion considers the general stress state by using the principal (difference) stresses, while the Dang Van 2 criterion determines a combination of shear and hydrostatic stress.

Solid mesh convergence checking

Previous results show that for some cases, application of Dang Van 2 criterion preserves a thin layer of non-eroded elements on the interface, and a kind of discontinuity for stress amplitude σ am DV 2 . This layer is particularly noticeable in Figures 4. 3 and 4.11. This singularity could be physically relevant, or the result of a numerical inaccuracy due to local mesh density. In order to check the validity of mesh near the interface, the height of elements has been divided by four (in the z direction). Then the results induced by both meshes are compared for one case (φ = 1 mm,V = 100 m.s -1 ). Figure 4.15 shows the comparison of the Dang Van 2 stress amplitude σ am DV 2 for both meshes. The result is obviously smoother using refined mesh, but the same tendency is observed near the interface. For more relevance, eroded zone are compared between both criteria in Figure 4.16 for N lim = 10 7 cycles. The fatigue function E DV 2 is also smoother with refined mesh, but the removed broken elements constitute globally the same erosion shape. Moreover, the layer is still present, and owns the same thickness. Finally, the comparison is set from a qualitative aspect by computing the eroded region area (since the model treats plane strains). Its value equals the sum of each element multiplied by it surface. As the affected region is composed of a homogeneous mesh, total area becomes simply the total number of elements times the common element surface. • The standard mesh is made of perfect squares with a a = 5 • 10 -6 m side. The number of eroded elements with this mesh is 29. The total erosion surface is then equals to 725 µm 2 .

• The refined mesh is composed of rectangular elements with width a, and height a/4. There are 122 eroded elements. Thus the total erosion area is 763 µm 2 .

• Finally, the error ratio between the two areas is almost 5%.

This ratio is low, and proves the convergence of mesh. Moreover the non-eroded layer seems to be physical, because it remains even for an element thickness four times lower.

Adding the coating layer

In this section, the coating layer described previously (in 4.1) is implemented on top of the solid domain. The comparison between the pure steel and the coated material is made in terms of transient evolution first, then in terms of lifetime by fatigue analysis. Both criteria are considered, but we restrict ourselves to the case (φ = 1 mm,V = 100 m.s -1 ).

Effects on transient evolution

In this section the effect of coating layer is studied in terms of time evolution of stresses without fatigue considerations. The evolution of hydrostatic and Von Mises stresses are given in Figures 4.17 = 4233 m.s -1 , which are definitely faster than raw material ones. This means that primary and secondary waves (and thus Rayleigh waves) propagate faster in the coating layer than in the raw material.

The different frames depicted in Figures 4.17 • Frame (a) : After impact, but before the P-waves reach the interface between coating and raw material.

• Frame (b) : After P-waves reach the end of coating, which generates reflected waves in the coating and refracted ones in the raw material.

• Frame (c) : After multiple reflections of waves inside the coating between the free surface and the interface with raw material domain.

• Frame (d) : After the (refracted) waves have travelled a more important part of the raw material domain.

On frame (a) we see clearly the time advance of the P-waves (and S-waves in a lesser extent) in the coating (right) ahead of the case of pure steel (left). This lag is amplified as the waves travel in the coating. The gap is therefore distinctly seen on frame (b) and mostly (c) in the raw material after refraction of the waves at the material interface. The Frame (d) shows that as predicted Rayleigh waves propagate faster in coating, which is particularly observable with hydrostatic stress in Figure 4.17. The Von Mises stress (Figure 4.18), which denotes the intensity of shear, seems to be subdued in the coating, because despite its advance in this model compared to the non-coated one, its value is reduced in the raw materiel at the same time, which is clearly shown on frame (d). The difference is less pronounced for the hydrostatic stress, what could be explained by the stiffness of the coating, which transmits more the compression effects than the shear ones because it undergoes less distortion than the raw material. of waves velocities between the two media, which allows a smaller part of elastic energy to cross the interface than the pure steel model. For Dang Van 2 criterion, we observe an increase of the stress σ am DV 2 in the coating near the interface with raw material, which does not appear for ASME criterion. This could be due to the reflection of waves on the interface and a higher shear, as the Dang Van 2 criterion dissociates the traction/compression and shear effects. Indeed, Figures 4.17 and 4.18 from previous section show that hydrostatic stress remains low (in absolute value) in this region, whereas the Von Mises shear, which is more uniform. The interface in this region is compressed on both sides (σ H < 0), but the thin region inside the coating is only sheared (σ H 0 and σ VM = 0).

Due to its higher resistance to fatigue, the coating has a higher number of cycles to failure (both criteria) than the raw material (Figures 4.20, 4.23 and 4.24), although the stress amplitudes have similar value. This leads to the erosion of raw material before the coating (Figures 4.25 and 4.26). But the coating diminishes significantly the erosion of raw material because it carries the most part of the load (stresses). The number of eroded elements for N lim = 3 • 10 9 cycles is equal to 3329 and 3424 according to ASME and Dang Van 2 criteria respectively without coating, and become 5 and 22 with the coating layer. Only the raw material is damaged, as the coating has a good resistance to fatigue.

Both criteria confirm the same influence of coating layer, since the eroded region locates at the same position, which is under the material interface, on the center of solid (x = 0). In terms of stresses the maximum value locates near the surface (in the coating) near the microjets formation (x R jet ), which justifies the use of coating layer. The first element of the raw material is broken at N ≈ 2.4 • 10 9 cycles for both criteria, which corresponds to L h = 16 000 h, and is almost equivalent to 2 years for a turbine operating 24/7, which is relevant.

Discussion about results and limitations

The fact that raw material breaks before coating could not be relevant, as the coating is usually removed before the raw material domain. But as previously mentioned, it has a better resistance to fatigue than the raw material (S-N diagram), so the reason could be the physical model itself, like the basic material interface condition used here, or the solid plain strain description which could underestimate the stresses inside coating. Moreover the coating manufacturing process usually generates residual stresses inside the coating layer, which are not taken into account in the present work. They could change the mean stress in the fatigue load and decrease the lifetime if traction effects are amplified. Even if raw material breaks first, a second computation with eroded elements in the raw region could lead to the failure of the coating.

The thin virgin elements layer found with Dang Van 2 criterion (cf. section 4.2.2) is not present using the ASME criterion due to their different formulation because each of them does not focus on the same physical quantities. This layer has a far higher number of cycles than the eroded region : N ≈ 10 9 vs. N ≈ 10 7 cycles according to fatigue analysis, but it would naturally break on the next physical impact, i.e. for N = N lim + 1, which is quite equivalent to N N lim , because N lim 1,i.e. one impact is negligible compared to the 10 7 previous ones.

As mentioned in Chapters 1 and 3, fatigue damage due to droplet impingement is caused by crack formation and propagation at mesoscopic scale. As the current model deals with a continuum isotropic homogeneous medium, when an element of the mesh is eroded it corresponds to the presence of a crack inside, but its direction and length are unknown, because the crack propagation itself is not modelled. This erosion method is effectively mesh-sensitive because the eroded zone is defined by the envelope of the phantom elements edges. The new FSI surface depends on the size and shape of elements.

These results are generated considering strong assumptions : there is no surface change during the N lim cycles and the surface is cleaned between each cycle, thus water is removed before the next impingement, and shall not influence the next one. Therefore, a small number of cycles N lim must be chosen for realistic model, because a different surface would certainly modify the physics of impact, such as stress concentration inside solid or higher fluid compression. Nevertheless, N lim must be high enough to agree with high-cycle fatigue theory.

Finally, the present method does not include the impact frequency, although this could be an influential parameter on the lifetime according to [OKA 07].

Conclusion on prediction of damage due to droplets impacts

In this chapter, the fatigue analysis method described in Chapter 3 is performed to simulate and analyse the damage. First, a parametric analysis for several droplet diameters and 127 impact velocities is investigated with both fatigue criteria in order to determine the influence of each parameter on the damage intensity. It appears that increasing them raises the damage intensity, and velocity has a far higher impact than diameter, mostly if both parameters are combined. For low velocities, the eroded regions locates clearly around the envelope of microjets area. Application of both criteria gives globally similar results, and the difference of results change depending on the cases. ASME criterion gives higher maximum stress amplitude so lower minimum number of cycles than Dang Van 2 one. Eroded area is larger according to ASME criterion for low velocity, whereas for high velocity, but erosion depth is almost the same for both criteria. Differences of results between them are probably due to their respective formulation and the physical quantities on which each one focuses. Then, the effects of a coating layer on the surface of solid are studied in terms of transient evolution and fatigue lifetime improvement. The coating being more stiff than the raw material it reduces the shear in this part, and the difference of waves velocities between both media generates reflection inside the coating and only a part of elastic energy is transmitted to the raw material. Concerning the fatigue analysis, the coating improves significantly the bucket lifetime. Indeed, its carries the same load as non-coated material, but has a better fatigue resistance than raw material. Moreover most part of stress amplitude is confined to the coating due to waves reflection. Adding the coating layer multiplies by 10 3 the lifetime, which corresponds to a more relevant value, since the coating is effectively used to increase the lifetime of Pelton buckets.

General conclusion & perspectives

The first chapter described the mechanisms responsible for erosion of metallic structure, as Pelton buckets by water droplet impingement. The first part deals with the fluidstructure interaction itself between the spherical water droplet and the plane solid target.

Concerning the fluid, at the beginning of impact, a shock wave travels from the contact zone to the opposite side. The free surface acts like a mirror and generates a reflected wave in the fluid. After reaching the top of droplet, the superposition of relaxation waves drops the pressure. The compression waves are present inside the solid in the same way, starting from the impact point, but are followed by shear waves too, and to a lesser extend by Rayleigh waves on the surface. The contact pressure can be firstly described by the water-hammer pressure, whose value can be analytically found, and denotes the pressure on the center of the contact at start of impingement. The maximum contact pressure, higher than waterhammer pressure, locates on the edge of contact area, and is found when microjets form. After that, its value diminishes to the stagnation pressure. These microjets are generated by the compression effects of the shock wave travelling inside the fluid. When the expanding contact area becomes subsonic, shock front overtakes it and compression effects eject the liquid in contact with solid. The velocity of microjets is far higher than the ambient sound speed of liquid. The reflected waves inside the droplet force local density to vanish, and this phenomenon could induce cavitation by phase change, but this behaviour will not be included to forthcoming simulations. For small droplet at high velocity, fluid behaviour is driven by inertial effects, thus viscosity and surface tension could be neglected for a relevant modelling.

The second part concerns the erosion mechanism itself, by experimental approach found in literature. The main test rigs for droplet impact erosion are presented, such as the Pulsating Jet Test Rig (PJET) or the Rotating Arm Apparatus, which is often used for rain impact erosion. Then, the macroscopic models and experimental studies on erosion are described. These models and studies give a general description of different phases observed : an incubation phase, followed by a linear erosion, and finally the last one where erosion rate decreases. Experimental investigations from literature analyse the influence of parameters like droplet velocity, its size or impact angle. It has been found that high velocity and size with normal impact increase damage. Next, the mesoscopic observation of the damage is presented through experimental works from literature. The mechanism is described by different steps. First, the first impacts erode grain boundaries and generate pits between grains. Then, microckacks appear at the bottom of these pits. After a lar-ger number of impacts, material is removed from surface when cracks join together. The erosion mechanism is finally fatigue cracking according to both macroscopic and mesocopic observations. This last information has been crucial as the basic mechanism for the development process of the erosion tool in Chapter 3.

In the second chapter, the numerical methods build for droplet impact simulation have been presented. The physics is based on continuum mechanics, thus the equations of solid and fluid mechanics have been presented, such as conservation laws and constitutive/state equations. Coupling conditions needed to be added to simulate correctly the interaction between the two sub-domains. Then, the numerical disretization methods have been described. The solid is discretized by the Finite Element Method and the fluid by the Smoothed Particle Hydrodynamics. To ensure a good compromise between accuracy and CPU time saving, a numerical convergence study is performed for the fluid computation. The quantity used to evaluate the validity of results was the pressure impulse, because it drives the response of the mechanical systems in fast transient dynamics. Results have shown that considering a two dimensional problem was relevant at this stage of the research. One only considers a symmetric model, and a truncation of the top on the droplet, depending on the simulation time and the impact velocity. The particle size has been investigated too, to find the best equilibrium. Next, a one-way computation has been build in two steps. A first simulation with the droplet impacting a rigid wall was performed, and the time distribution of contact pressure on wall was recorded. The second simulation involves the structure, and the previous pressure is entered as a load on the interface edges of the elements. This procedure allows understanding the general behaviour the whole system, but is complex to implement, and not adaptable to modifications of interface shape. On this basis, considering a two-way coupling should be preferred. This last simulation is a stronger model for analysing the structure behaviour. The stresses involved are not high enough to cause a sudden fracture of the material. Nevertheless, the regions near the contact are loaded by alternating traction and compression, which is a predominant factor for a fatigue cracking. This simulation allowed corroborating the fatigue cracking as damage mechanism for the current droplet impact phenomenon.

The third chapter deals with fatigue developments. First the fatigue phenomenon is presented in a global point of view, and the main notions and quantities significations are explained. Then the different classifications of fatigue as function of stresses and type of loads are distinguished. The 2-way coupled transient computation of Chapter 2 deals with high-cycle fatigue with non-proportional load. These properties allow choosing relevant fatigue criteria : the general ASME criterion, which is built from general considerations but is easy to implement, and the Dang Van 2 criterion. These criteria provide a local condition of erosion with the value of fatigue function. This analysis is useful for the current erosion simulation. Indeed the procedure consists of a post-processing computation after the FSI simulation. This operation performs a calculus on every element of the solid mesh using the previous fatigue criteria, its result allowing finding if each element is eroded or not for a given number of load cycles. And this condition of erosion is determined by the fatigue criteria. Eroded elements must be converted to "ghost", i.e they are not considered in next simulations.

In the last chapter, the fatigue analysis method described in previous chapter has been performed to simulate and analyse the damage. First, a parametric analysis for several droplet diameters and impact velocities is investigated with both fatigue criteria in order to determine the influence of each parameter on the damage intensity. It appears that increasing them raises the damage intensity, and velocity has a far higher impact than diameter, mostly if both parameters are combined. For low velocities, the eroded regions locates clearly around the envelope of microjets area. Application of both criteria gives globally similar results, and the difference of results change depending on the cases. ASME criterion gives higher maximum stress amplitude so lower minimum number of cycles than Dang Van 2 one. Eroded area is larger according to ASME criterion for low velocity, whereas for high velocity, but erosion depth is almost the same for both criteria. Differences of results between them are probably due to their respective formulation and the physical quantities on which each one focuses. Then, the effects of a coating layer on the surface of solid have been studied in terms of transient evolution and fatigue lifetime improvement. The coating being more stiff than the raw material it reduces the shear in this part, and the difference of waves velocities between both media generates reflection inside the coating and only a part of elastic energy is transmitted to the raw material. Concerning the fatigue analysis, the coating improves significantly the bucket lifetime. Indeed, its carries the same load as non-coated material, but has a better fatigue resistance than raw material. Moreover most part of stress amplitude is confined to the coating due to waves reflection. Adding the coating layer multiplies by 10 3 the lifetime, which corresponds to a more relevant value, since the coating is effectively used to increase the lifetime of Pelton buckets.

The method used and the tool built in this work have the advantage of possibility working for other types of damage sources, such as other PrEDHyMa projects : the hydroabrasive and the cavitation erosions, since they could be based on fatigue too.

The outlook for further research emerging from the current work is vast. As presented in Chapter 3, the "two-steps" simulations with loop of a second transient computation and fatigue analysis after a first one could be more relevant than the current models. Indeed, the change of interface shape and the presence of initial damage due to the preceding impacts should definitely modify the droplet behaviour and induce a different stress intensity and distribution.

Then the coating layer usually contains residual stresses as mentioned in Chapter 4, which have been neglected in this study. These stresses should change the fatigue resistance of the coating and provide more realistic results.

The current model considering an isotropic, homogeneous and continuous solid, this modelling lacks accuracy, but constitutes a good basis for the erosion simulation. As the damage is caused by intergranular fractures (Chapter 1), a more elaborate simulation could consider a mesoscopic approach, by modelling the behaviour of grains and joints.

Next, the solid computation deals with plane strains, which may not be realistic for the fatigue analysis particularly. Indeed, the stress state could be underestimate as well as overestimated. But a full 3D computation with both fluid and solid domains was limited anyway by the computations equipments of the laboratory.

The strong assumptions considered in this work on the droplet size and velocity should be replaced by a more realistic approach. Indeed, the fact that the same droplet impacts the bucket at the same velocity, on the same location, and at the same frequency must not be relevant. These features could be determined by statistical calculations to simulate the desired damage.

This fatigue damage model can also be applied directly to cavitation erosion as well as sand or solid erosion mechanisms simulation.

Finally, as the release waves inside the droplet induce a pressure drop, cavitation could form and become an influential phenomenon in droplet impact erosion. For this reason a phase-change model for the fluid could be inserted in the fluid behaviour to combine the action of impact and cavitation, and investigate the effects of each damage source. Au moment de l'impact, une onde de compression se progage dans la goutte à vitesse c f ≥ c 0 f en partant du point de contact vers l'autre extrémité (Figure 27). Elle est de plus réfléchie élastiquement le long des surfaces libres par effet miroir. Au même moment, la zone d'impact évolue et croît en taille, ce qui a pour conséquence le changement de forme du front d'onde. Il est bien-entendu possible de construire géométriquement ce front à l'aide du principe de Huygens-Fresnel [FIE 99]. Le front est ainsi l'enveloppe des fronts d'ondelettes générées par les bords successifs du contact. Le volume contenu à l'intérieur du front est en forte compression, et sa masse volumique est supérieure à la valeur initiale ρ f > ρ 0 f . Concernant le solide, l'impact de la goutte provoque aussi une onde de compression longitudinale (onde P) dans sa masse, de vitesse c L , mais aussi une onde de cisaillement transverse (onde S) de vitesse c T < c L . Ces deux ondes sont de type sphérique, mais on observe aussi la propagation d'une onde de surface ou onde de Rayleigh.

Explication physique de l'érosion par choc de goutte 0.1.1.2 Pression de contact et microjets Lors de l'impact d'un corps fluide sur un solide, on observe le bien connu effet "coup de bélier". Ce phénomène génère une pression p wh appelée pression coup de bélier donnant la valeur de la pression de contact p(x,t) au centre du contact et au moment exact de l'impact :

p(x = 0,t = 0) = p wh (4)
Pour un mur solide rigide, [FIE 12, HEY 69, KEN 00] proposent l'expression :

p wh = ρ 0 f c f V (5)
La vitesse c f peut être approximée par la relation linéaire c f = c 0 f + k f V pour une vitesse d'impact inférieure à 1000 m.s -1 . La constante k f est liée au fluide et vaut environ 2 pour l'eau. Si l'on observe la distribution temporelle de la pression de contact, on constate qu'à chaque instant, le maximum se situe toujours au bord de la zone de contact. La valeur maximale apparait au moment où le front d'onde dans la goutte dépasse la zone de contact. A ce moment, les effets de compression avec la cible provoquent l'éjection latérale du fluide le long de la zone de contact en microjets (Figure 28). Ainsi le pic maximal en pression p max est obtenu à l'instant et l'endroit d'apparition des microjets :

p max = p(x = R jet ,t = T jet ) (6)

Caractéristiques du modèle numérique

La goutte d'eau est un disque de rayon R = 0.5 mm impactant normalement le solide à vitesse V = 100 m.s -1 , de masse volumique initiale ρ 0 f = 1000 kg.m -3 , vitesse du son initiale c 0 f = 1500 m.s -1 sans viscosité ni tension superficielle. Le caclcul est réalisé en bidimensionnel et bi-directionnel (calcul plan) car la méthode SPH ne connait actuellement pas de description axisymétrique. L'absence d'effet convectif permet de choisir la loi de Tait isotherme barotrope comme loi d'état du fluide, donnée par l'équation (7), où γ f = 7 représente le coefficient polytropique du fluide :

p f = ρ 0 f (c 0 f ) 2 γ f   ρ f ρ 0 f γ f -1   (7)
Le domaine solide est alors un rectangle d'acier X3CrNiMo13-4 homogène isotrope considéré élastoplastique à écrouissage isotrope, de module d'Young E 0 = 200 GPa, module tangent E t = 20 GPa, coefficient de Poisson ν = 0.288, limite élastique initiale σ 0 Y = 560 MPa et de masse volumique ρ 0 s = 7700 kg.m -3 . Les dimensions du rectangle (H, 2H) sont sufisamment élevées pour éviter la réflexion des ondes P (les plus rapides), et la théorie des déformations planes est considérée.

Le tout est représenté sur la Figure 30. L'intégration temporelle est réalisée par un schéma différent dans chaque domaine, mais explicite dans les deux cas : le domaine fluide est discrétisé par le schéma de Runge-Kutta 2 (ou Midpoint), et le domaine solide par celui de Newmark.

Equations régissant le comportement des sous-domaines

Les quantités exprimées simultanément dans les deux sous-domaines sont différenciées par les indices suivants : • f pour le fluide et • s pour le solide. Le déplacement est noté u, la vitesse v, la contrainte de Cauchy σ, la déformation ε, la masse volumique ρ, la pression p, la matrice identité 1 et les forces volumiques à distance b.

Le domaine fluide est décrit par les équations d'Euler données ici en forme Eulérienne et l'équation de Tait :

                       ∂ρ f ∂t + ∇ x • (ρ f v f ) = 0 ∂(ρ f v f ) ∂t + ∇ x • (ρ f v f ⊗ v f ) + ∇ x • (p f 1) = ρ f b p f = ρ 0 f (c 0 f ) 2 γ f   ρ f ρ 0 f γ f -1   (8)
Pour le solide, l'équation classique de conservation de la quantité de mouvement est utilisée (ici exprimée en configuration Eulérienne). Le matériau étant considéré élastoplastique, les déformations ε s sont partitionnées en une déformation purement élastique ε e s et une purement plastique ε p s . La partie plastique des déformations n'est pas détaillée dans ce résumé. L'ensemble est décrit dans l'équation suivante :

                     ∇ x • σ s + ρ s b = ρ s ∂ 2 u s ∂t 2 ε s = 1 2 ∇ x u s + ∇ x u s + ∇ x u s • ∇ x u s ε s = ε e s + ε p s σ s = λtr (ε e s ) 1 + 2µε e s (9)

Convergence et impact sur mur rigide

Des calculs pour étudier la convergence du modèle en comparant le chargement induit par la goutte sur un mur rigide ont permis de déterminer les caractéristiques optimales de la modélisation pour gagner en rapidité, sans perdre en précision ni stabilité. Il s'agit de la taille des particules, la hauteur de troncature du domaine fluide (sur (O, z)),la symétrie plane (normale à x) et la possibilité de calcul 2-D. La quantité physique considérée pour la comparaison est l'impulsion, définie comme l'intégrale temporelle de la force pendant une durée choisie. En effet, pour un oscillateur de masse et raideur unitaires, sans amortissement (négligeable en dynamique rapide) excité par une force constante F 0 durant un temps court T F , la réponse en déplacement z + (t) et en vitesse ż+ (t) après relachement est la suivante :

Pour T F 1,          z + (t) = F 0 T F sin(t) - T F 2 cos(t) ż+ (t) = F 0 T F cos(t) + T F 2 sin(t) (10)
Or le produit F 0 T F est simplement l'impulsion de F 0 durant le temps T F . On voit donc dans l'équation [START_REF] San | [END_REF] que le déplacement et la vitesse sont pilotés par l'impulsion en mécanique des chocs.

La comparaison de l'impulsion permet de montrer notamment la pertinence de l'utilisation d'une symétrie, de la possibilité d'une troncature (hauteur dépendant du temps de calcul et de la vitesse d'impact), du passage à un modèle 2-D, et de déterminer la taille maximale des particules SPH.

La simulation d'impact sur un mur droit rigide est utile pour la compréhension du type de chargement induit par la goutte sur le solide. La Figure 31 montre la répartition spatiale de la pression en différents instants, et une pression de contact maximale p max = 170 MPa, équivalente à ≈ 10p wh . Ce pic de pression est localisé à x/R 0.18, soit dans la même zone que cas étudié par [HAL 02b]. L'impulsion par unité de surface donne une meilleur idée de l'intensité du chargement. La valeur maximale est au centre de l'impact (x = 0) et des fluctuations peuvent être observées en x/R = 0.18, une instabilité locale étant provoquée par l'apparition des microjets.

Simulation par couplage faible

Un premier calcul en couplage faible a été réalisé : la distribution spatiale de pression est relevée au cours du temps sur le mur rigide (simulation sur le domaine fluide uniquement), puis dans un second temps, cette pression est utilisée comme chargement sur le solide (simulation sur le domaine solide uniquement). Cette méthode est limitée pour plusieurs raisons : une grande précision spatiale et temporelle est demandée dans l'acquisition du chargement, les valeurs proviennent du calcul sur mur rigide, ce qui ne correspond pas au comportement du solide, et cette technique nécessite deux calculs à la suite. Pour toutes ces raisons, un calcul couplé fort est préféré.

Simulation par couplage fort

L'analyse des contraintes dans le solide permet de mettre en lumière les différents phénomènes se produisant dans le milieu. La contrainte hydrostatique σ H (x,t) permet d'observer la propagation des onde P et de connaître l'état de traction/compression du matériau en tout point x et temps t. La contrainte de Von Mises σ VM (x,t) concerne les ondes S, et donne l'intensité du cisaillement. La Figure 32 délimités, ceci étant dû à la zone de contact qui évolue dans le temps, provoquant des sollicitations combinées, et de plus à grande vitesse, ce qui a tendance à estomper la frontière entre les fronts.

L'évolution temporelle de ces deux quantités est donnée sur la Figure 33. L'analyse de la triaxiatlité des contriantes ζ montre que le matériau subit la combinaison des deux sollicitations.

Les valeurs des contraintes sont faibles comparées à la limite d'élasticité σ 0 Y , ce qui laisse penser que le mécanisme d'érosion est la fatigue à grand nombre de cycles, ceci permettant de choisir le futur critère de fatigue. La contrainte de Von Mises signée σ sign VM dont la valeur indique l'intensité du cisaillement et le signe le sens de l'effort sphérique (traction/compression) est un outil pertinant dans la détermination des zones sollicitées par la fatigue, car les fissures s'amorcent dans les zones cisaillées, et progressent lors d'un état de traction (la compression les bloquant). On étudie aussi l'évolution des sontraintes principales (σ I , σ II , σ III ) dont la plus grande valeur indique l'état maximal de contrainte en traction/compression, et les valeurs sont du même ordre de grandeur que les contraintes précédentes. Ensuite, l'observation des rapports de contraintes principales r i, j et de l'angle de rotation du repère principal α p indique un chargement non proportionnel (non radial), ce qui est déterminant dans le choix du critère de fatigue.

Procédure d'érosion

La simulation d'érosion est réalisée par post-traitement à la suite du calcul transitoire couplé. Pour chaque élément (un seul point de Gauss), on calcul une fonctin de fatigue E ∈ [0, 1] dont la valeur détermine l'état de l'élément. Lorsque E = 1, l'élément est considéré comme érodé et devient "fantôme", par contre tant que E < 1, l'élément reste dans le maillage, la valeur E = 0 correspondant à un élément vierge. La fonction de fatigue est en fait le rapport entre une amplitude de contrainte σ am et une contrainte limite σ lim provenant de la courbe de fatigue S-N (ou de Wöhler) du matériau. En effet, la fonction de fatigue est donnée pour un nombre de cycles limite N lim choisi par l'utilisateur, dont la contrainte correspondante par la courbe de fatigue est σ lim = f SN (N lim ) :

E(N lim ) = σ am σ lim (11)
L'amplitude de contrainte calculée peut être aussi reportée pour donner le nombre de cycles à la rupture N = f -1 SN (σ am ). La condition E(N lim ) = 1 correspond donc à N ≤ N lim autrement dit à σ am ≥ σ lim . Après le premier calcul couplé avec érosion, la nouvelle interface FSI est redéfinie, et un autre calcul couplé avec post-traitement peut être réalisé. La procédure est résumée sur la Figure 34. L'amplitude de contrainte σ am est calculée à partir de critères de fatigue décrits ciaprés.

Critères de fatigues utilisés

On utilise dans ce travail deux critères différents : un premier critère général proposé par l'American Society of Mechanical Engineers (ASME). Le second est le critère à plan critique de Dang Van 2 (DV2). Les critères à plan critique considèrent que la contrainte d'amplitude maximale se situe sur un plan de l'espace particulier, le plan critique. Ce critère dissocie les effets du cisaillement qui détermine la condition d'amorc ¸age des fissures d'une part, et la contrainte hydrostatique qui pilote la condition de leur propagation. Les termes relatifs à chaque critère sont notés avec les indices suivants : • AS pour le critère ASME et • DV 2 pour le critère de Dang Van L'amplitude de contrainte σ am AS est calculée à partir des contraintes principales, car elles donnent l'état de contrainte maximal dans le matériau. L'amplitude est la moitié de l'étendue de contrainte ∆σ, calculée comme le plus grand écart de contraintes au cours du cycle. On considère deux temps physiques t m et t n de containtes de Cauchy respectives σ (m) et σ (n) . On définit alors le tenseur des différences ∆σ (m,n) pour le couple (t m ,t n ) comme :

∀ (m, n) ≤ n T , m = n ∆σ (m,n) = σ (m) -σ (n) (13)
où n T est le nombre de pas de temps considérés dans le calcul de la contrainte équivalente. La diagonalisation de ∆σ (m,n) conduit aux différences de contraintes principales ∆σ La fonction de fatigue de ce critère peut être écrite comme suit :

E DV 2 = max t∈[0,T c ] ¨τa (t) + α DV σ H (t) β DV « (15)
où σ H (t) est la contrainte hydrostatique, τ a (t) le cisaillement alterné, et α DV et β DV sont deux constantes liées au matériau telles que :

α DV = 3 τ -1 σ -1 - 1 2
β DV = τ -1 (16)
où σ -1 and τ -1 sont les limites d'endurance en traction et torsion symétriques alternées respectivement. Les mesures permettant l'obtention de τ -1 étant complexe à réaliser, la valeur approximative τ -1 = σ -1 / √ 3 est utilisée, car elle se montre adaptée à l'acier considéré. De plus ce choix satisfait à la condition du critère α DV > 0. On note que le paramètre α DV ne dépend plus du matériau.

Pour pouvoir mettre la fonction de fatigue de l'équation (15) sous la forme standard de l'équation (12), on définit l'amplitude de contrainte comme suit, avec σ lim = σ -1 : L'étude paramétrique concerne les valeurs de vitesse d'impact V = {100 m.s -1 , 200 m.s -1 } et les diamètres φ = {0.5 mm, 1 mm, 2 mm}. Pour chacun des deux critères, l'amplitude de contrainte σ am , le nombre de cycles à la rupture N, la fonction de fatigue E(N lim ) pour un nombre de cycle N lim et les caractéristiques de la zone érodée sont comparés entre les différents cas (φ i ,V j ). Dans touls les cas, les contraintes sont toujours en dessous de la limite élastique du matériau σ 0 Y , ce qui justifie l'hypothèse de fatigue à grand mombre de cycles. On voit aussi que le critère ASME est plus sévère que celui de DV2 pour la contrainte maximale sur tout le maillage, et cet écart s'accentue avec l'augmentation de V . Le cas engendrant la plus grande contrainte pour les deux critères est (φ = 2 mm,V = 200 m.s -1 ), ce qui parait logique compte tenu de la valeur de l'énergie cinétique de la goutte avant impact. La vitesse est un paramètre plus influent que le diamètre de la goutte dans les valeurs de contraintes, et donc le nombre de cycles à la rupture. Dans le cas (1 mm, 100 m.s -1 ) la zone avec les plus hautes contraintes se situe à l'extérieur de la formation théorique des microjets.

σ am DV 2 = √ 3
Bien que σ am < σ 0 Y dans tous les cas, l'étendue des nombres de cycles à la rupture va de N = 10 4 pour le cas (2 mm, 200 m.s -1 ) à 10 6-7 pour (0.5 mm, 100 m.s -1 ), ce dernier cas correpondant presqque à de la fatigue gigacyclique.

La condition d'érosion choisie est N lim = 5 • 10 7 cycles, ce qui permet d'obtenir une zone d'érosion pour tous les cas (φ i ,V j ). L'intensité du dommage est caracérisée à l'aide des paramètres suivants : le nombre total d'éléments érodés, la profondeur d'érosion, la forme de la zone érodée et sa position.

Une fois encore, la vitesse d'impact se montre le paramètre le plus influent pour tous ces indicateurs, mais le diamètre contribue tout de même à l'érosion. Le nombre d'éléments peut être multiplié par 100, voire par 1000 pour une vitesse double. La profondeur d'érosion se voit multipliée par 5 pour φ = 2 mm.

Pour les faibles vitesses, la zone de dommage est plutôt localisée sous le contact, alors qu'elle s'étend loin près de la surface pour les vitesses plus élevées. Les deux critères montrent un dommage situé sous les zones de formation des microjets, en forme de couronne, visible sur la Figure 35.

Ajout de la couche de coating

On évalue maintenant l'effet d'une couche de coating en carbure de tungstène sur le matériau de base pour plus de réalisme. Le coating est modélisé comme un milieu continu, homogène et isotrope d'épaisseur 300 µm distinct du domaine constitué du matériqu de base, avec pour masse volumique ρ 0,coat s = 13500 kg.m -3 , module d'Young E coat 0 = 600 GPa, coefficient de Poisson ν coat = 0.24 et limite d'élasticité σ 0,coat Y = 1.11 GPa. Un seul cas est considéré pour la comparaison : (φ = 1 mm,V = 100 m.s -1 ).

D'abord, les effets sur l'évolution temporelle du système sont traités. Pour cela on observe les contrainte hydrostatique σ H et de Von Mises σ VM . Le coating possède des vitesses de propagation des ondes plus élevées que le matériau de base. On voit donc un décalage des fronts d'ondes entre les deux modèles. La contrainte de Von Mises se retrouve réduite lors de la présence de coating (Figure 36), car le coating est bien plus raide (module d'Young) que le matériau de base, et transmet ainsi moins le cisaillement dans la matière.

Les effets sur le dommage en fatigue sont aussi abordés, et on compare de nouveau l'amplitude de contrainte σ am , le nombre de cycles à la rupture N, et la fonction de fatigue E(N lim ), cette fois pour N lim = 3 • 10 9 cycles, car le coating améliore significativement la résistance à la fatigue de la structure (Figure 37). Le coating connaît les mêmes contrain- tes que le modèle simple, mais elles se retrouvent alors réduites dans le matériau de base. Elles sont divisées par quatre, et le nombre de cycles à la rupture est multiplié par ≈ 10 3 , les deux critères donnant des résultats semblables. Ce changement de contrainte peut s'expliquer par la différence de vitesses du son entre le coating et le matériau de base, et donc qu'une partie de l'énergie élastique est transmise par réfraction, l'autre partie étant réfléchie. L'application du critère de DV2 provoque l'apparition d'une concentration de contraintes dans le coating près de la zone de transition des deux matériaux, ce qui n'apparait pas dans l'utilisation du critère ASME. Cela peut être engendré par un fort cisaillement à l'interface, comme le critère de DV2 dissocie le cisaillement de la contrainte hydrostatique, contrairement au critère ASME. En effet, on observe une réduction de la contrainte hydrostatique à l'interface, alors que la contrainte de Von Mises ne montre pas de gradient local. La présence de coating reporte la rupture de l'acier à N ≈ 2.4 • 10 9 cycles, soit après environ 2 ans de fonctionnement 24/7 de la turbine pour les deux critères, ce qui reste pertinent, et montre bien l'intérêt de l'application d'une couche de coating. Les simulations numériques par intéraction fluide-structure qui suivent permettent de simuler le dommage par choc de goutte. Le domaine solide est discrésisé par la MEF, et le fluide par la méthode SPH. Après une étude de convergence permettant de réduire la taille du problème, un calcul transitoire couplé a confirmé l'hypothèse de dommage par fatigue en analysant les états de contraintes.

Conclusion et perspectives

L'érosion en soit est traduite à l'aide de critères de fatigue qui fournissent une condition sur chaque élément du maillage solide pour l'éliminer, ce qui revient à représenter une perte de matière. Ce calcul est effectué par post-traitement à la suite du calcul transitoire, où un nombre de cycles virtuels est imposé. Les critères de fatigue retenus sont au nombre de deux : le critère général ASME pour sa simplicité d'implémentation, et le critère de Dang Van 2, adapté à la fatigue à grand nombre de cycles.

Une étude paramétrique a ensuite été réalisée pour différents diamètres et vitesses de gouttes. Il apparaît que la vitesse a une plus grande influence que le diamètre en termes de contraintes équivalentes, qui restent sous la limite élastique. Pour les vitesses faibles, on constate particulièrement les effets des microjets car la zone d'érosion se situe autour de leur formation.

Pour plus de réalisme, une couche de revêtement (coating) a été incluse dans le modèle, et une comparaison avec le modèle simple a été effectuée. Le coating par sa grande raideur ralentit les ondes P et S, et une réduction du cisaillenent a été détéctée en dessous de sa frontière, due à une transmission moindre des distorsions. Comme supposé, le coating a bien amélioré la durée de vie du matériau, multipliée par 1000. La méthode présentée a l'avantage de pouvoir être utilisée pour les autres types d'érosion du projet PrEDHyMa. Pour améliorer cette démarche, un bouclage pourrait être effectué après le premier calcul d'érosion. Ensuite, le modèle utilisé pourrait être amélioré, car il considère des hypothèses fortes, comme l'absence de contraintes résiduelles dans le coating, l'aspect continu, homogène et isotrope du matériau, le calcul en déformations planes, la similitude des gouttes impactantes, ou encore l'absence de cavitation dans le domaine fluide.

  Cette thèse débutée en octobre 2014 a été réalisée dans le cadre du projet européen FP7 / 2007-2013 sous la convention de subvention 608393, nommé PrEDHyMa. Elle fut le fruit d'une collaboration entre le Laboratoire de Mécanique des Contacts et des Structures (LaMCoS) de l'Institut National de Sciences Appliquée (INSA) de Lyon, le Laboratoire de Mécanique des Fluides et d'Acoustique (LMFA) de l' École Centrale de Lyon (ECL) et la société ANDRITZ Hydro. Je tiens donc à remercier en premier lieu l'ensemble des partenaires m'ayant fourni un cadre de travail propice à la réalisation de ce travail durant ces trois années. Je remercie MM. Vincent Faucher et Habibou Maitournam pour avoir accepté d'être rapporteurs de mon travail, Mme. Catherine Gardin pour avoir accepté de présider ce jury, ainsi que les examinateurs : M. Jean-Christophe Marongiu et Mme. Emmanuelle Sallé.

Introduction1

  Physical phenomenon of droplet impact erosion 1.1 Droplet impact . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.1 Waves propagation . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.2 Contact pressure . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.3 Jetting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.4 Cavitation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Erosion mechanism : experimental approach . . . . . . . . . . . . . . . .

  injector vertical-axis Pelton runner . . . . . . . . . . . . . . . . . . 7 2 Erosion damage visible on a Pelton turbine bucket at the end of its functional lifetime . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.1 Shock front and highly compressed volume on early time of droplet impact 12 1.2 Geometrical construction of wave front with Huygens-Fresnel principle . 13 1.3 Time evolution of density during the droplet impact . . . . . . . . . . . . 14 1.4 Different types of waves involved when a water droplet impacts an elastic, isotropic, and homogeneous solid body . . . . . . . . . . . . . . . . . . . 15 1.5 Water-hammer pressure vs. impact velocity for water droplet on rigid target 16 1.6 Ratio between water hammer edge and droplet radius vs. impact velocity for water droplet on rigid target . . . . . . . . . . . . . . . . . . . . . . . 17 1.7 Spatial distribution of contact pressure for several times after impact (R = 0.1 mm, V = 500 m.s -1 ) . . . . . . . . . . . . . . . . . . . . . . . . . . 18 1.8 Map of high pressure duration as a function of droplet radius and impact velocity for water droplet on rigid target . . . . . . . . . . . . . . . . . . 20 1.9 Birth of lateral microjets . . . . . . . . . . . . . . . . . . . . . . . . . . 21 1.10 Radial velocity and density distribution on the contact zone when jetting . 21 1.11 Steps of near-side cavitation formation . . . . . . . . . . . . . . . . . . . 22 1.12 Impact of 10 mm diameter water drop by a metal slider at a velocity of 110 m.s -1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 1.13 Principles scheme of PJET . . . . . . . . . . . . . . . . . . . . . . . . . 24 1.14 Example of a small, relatively low-speed, rotating disk-and-jet repetitive impact apparatus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 1.15 Example of a large, high-speed, rotating arm-and-spray distributed impact apparatus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 1.16 Typical erosion curve . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 1.17Normalized volume loss of X20Cr13 as a function of droplet impact angle with impact speed V = 488 m.s -1 . . . . . . . . . . . . . . . . . . . . . 27 1.18 Erosion behaviour of X20Cr13 at different droplet sizes tested at the impact velocity of 488 m.s -1 . . . . . . . . . . . . . . . . . . . . . . . . . 28 1.19 Erosion behaviour of 12%Cr steel at different impact velocities for a droplet size φ = 0.2 mm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 v 1.20 Grain ruptures on the initial surface of H300 after 50'000 impacts at impact velocity V = 225 m.s -1 . . . . . . . . . . . . . . . . . . . . . . . . 1.21 Fatigue defects observed on H300 after 5 million and 10 million impacts, at impact velocity V = 225 m.s -1 . . . . . . . . . . . . . . . . . . . . . 1.22 Cross section of eroded stainless steel at different magnification, at impact velocity V = 225 m.s -1 . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.23 Intergranular damage observed on Ti-6Al-4V, at impact velocity V = 350 m.s -1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.24 Born of triple junctions an microvoids on Ti-6Al-4V, at impact velocity V = 350 m.s -1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.25 Grain ejection, striation inside crater and slip bands on Ti-6Al-4V, at impact velocity V = 350 m.s -1 . . . . . . . . . . . . . . . . . . . . . . . . 1.26 Stress wave propagation observed on γ-TiAl surface . . . . . . . . . . . . 1.27 Intergranular and intragranular cracking in γ-TiAl . . . . . . . . . . . . . 1.28 Triple split resulting from grain pulled-out in γ-TiAl . . . . . . . . . . . . 1.29 The two modes of intragranular fractures in γ-TiAl . . . . . . . . . . . . 1.30 Areas sensitive to cracking in Pelton buckets . . . . . . . . . . . . . . . . 1.31 HVOF coating on Pelton buckets . . . . . . . . . . . . . . . . . . . . . . 1.32 Schematic diagram of the high velocity oxy-fuel spray process (HVOF) . 1.33 WC 10Co 4Cr coating materials . . . . . . . . . . . . . . . . . . . . . . 1.34 Section of a coated bucket . . . . . . . . . . . . . . . . . . . . . . . . . 1.35 Eroded surfaces of Pelton runner without coating and with SHX TM coating after operation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1 The normal vectors for each sub-domain at the fluid-structure interface . . 2.2 The solid domain Ω s discretized by a Lagrangian mesh under two types of boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Kernel approximation with the smoothing kernel function . . . . . . . . . 2.4 Representation and dimensions of the numerical model . . . . . . . . . . 2.5 Dynamic system and time evolution of load force . . . . . . . . . . . . . 2.6 Distribution of probes in 2-D model . . . . . . . . . . . . . . . . . . . . 2.7 Distribution of probes in 3-D model . . . . . . . . . . . . . . . . . . . . 2.8 Impulse distribution of 3-D computation . . . . . . . . . . . . . . . . . . 2.9 Comparison of impulse distribution between 2-D and 3-D models . . . . 2.10 Comparison of impulse distribution between full and symmetric models . 2.11 Different models for truncation . . . . . . . . . . . . . . . . . . . . . . . 2.12 Comparison of impulse distribution between full and truncated models . . 2.13 Comparison of impulse distribution for different particle sizes . . . . . . 2.14 Contact pressure on rigid wall . . . . . . . . . . . . . . . . . . . . . . . 2.15 Pressure field inside the droplet at t = 150 ns after impact with ∆x = 6 • 10 -7 m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.16 Hydrostatic stress inside the solid at t = 600 ns after impact . . . . . . . . vi 2.17 Pressure inside fluid, hydrostatic and Von Mises stresses inside solid for several times after impact . . . . . . . . . . . . . . . . . . . . . . . . . . 80 2.18 Time evolution of hydrostatic and Von Mises stresses at x = 0.1425 mm, z = -2.5 µm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 2.19 Time evolution of stress triaxiality at x = 0.1425 mm, z = -2.5 µm . . . . 83 2.20 Signed Von Mises Stress inside solid and pressure inside fluid at jetting time 84 2.21 Time evolution of principal stresses at x = 0.1425 mm, z = -2.5 µm . . . 85 2.22 Time evolution of principal stresses ratios at x = 0.1425 mm, z = -2.5 µm 86 2.23 Time evolution of principal angle at x = 0.1425 mm, z = -2.5 µm . . . . 87 3.1 Example of sinusoidal load cycles . . . . . . . . . . . . . . . . . . . . . 92 3.2 S-N diagram and different fatigue domains . . . . . . . . . . . . . . . . . 94 3.3 Decomposition of stress vector in basis linked to physical plane . . . . . . 98 3.4 Definition of terms related to normal stress . . . . . . . . . . . . . . . . . 98 3.5 Definition of terms related to tangential stress . . . . . . . . . . . . . . . 100 3.6 Endurance domain and two typical load paths . . . . . . . . . . . . . . . 102 3.7 Procedure for erosion simulation . . . . . . . . . . . . . . . . . . . . . . 105 3.8 S-N diagram with upper stress and lower cycle limits . . . . . . . . . . . 105 3.9 Procedure for calculation of stress amplitude for the ASME criterion . . . 107 3.10 Procedure for calculation of stress amplitude for the Dang Van 2 criterion 108 3.11 Example of mesh with eroded elements and modification of FSI interface 109 4.1 Representation and dimensions of numerical model involving coating layer 113 4.2 Raw and coating materials behaviour curves : Traction curve and S-N diagram . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 4.3 Stress amplitude inside solid according to ASME criterion on left and

  2 criterion between pure basic material and coated one at t = 460 ns after impact for φ = 1 mm and V = 100 m.s -1 . . . . . . . . . . . . . . . . . . . . . 4.25 Comparison of fatigue function for N lim = 3 • 10 9 cycles according to ASME criterion between pure basic material and coated one at t = 460 ns after impact for φ = 1 mm and V = 100 m.s -1 . . . . . . . . . . . . . . . viii 4.26 Comparison of fatigue function for N lim = 3 • 10 9 cycles according to Dang Van 2 criterion between pure basic material and coated one at t = 460 ns after impact for φ = 1 mm and V = 100 m.s -1 . . . . . . . . . . . 138 Front d'onde et volume comprimé en début d'impact . . . . . . . . . . . 176 Formation des microjets latéraux . . . . . . . . . . . . . . . . . . . . . . 177 Dommage par fatigue observé sur de l'acier H300 après 5 et 10 millions d'impacts à vitesse V = 225 m.s -1 . . . . . . . . . . . . . . . . . . . . . 179 Représentation et dimensions du modèle numérique . . . . . . . . . . . . 180 Pression de contact sur le mur rigide . . . . . . . . . . . . . . . . . . . . 183 Pression dans le fluide, contraintes hydrostatique et de Mon Mises dans le solide en différents instants après impact . . . . . . . . . . . . . . . . 184 Évolution temporelle des contraintes hydrostatique σ H et de Von Mises σ VM en x = 0.1425 mm, z = -2.5 µm . . . . . . . . . . . . . . . . . . . 185 Procédure de la simulation d'érosion . . . . . . . . . . . . . . . . . . . . 186 Fonction de fatigue et éléments non-érodés pour N lim = 5 • 10 7 cycles d'après chaque critère, au temps t = 460 ns après impact, pour le cas (φ = 1 mm,V = 100 m.s -1 ) . . . . . . . . . . . . . . . . . . . . . . . . 190 Comparaison de la contrainte de Von Mises entre le modèle en acier pur et celui avec coating en différents instants pour le cas (φ = 1 mm,V = 100 m.s -1 ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 191 Comparaison de la fonction de fatigue pour N lim = 3 • 10 9 cycles d'après le critère de Dang Van 2 entre le modèle en acier pur et celui avec coating à t = 460 ns après impact pour φ = 1 mm et V = 100 m.s -1 . . . . . . . . 192 ix List of Tables 2.1 Material data for numerical simulations . . . . . . . . . . . . . . . . . . 2.2 CPU time for different 2-D fluid models . . . . . . . . . . . . . . . . . . 4.1 Material data for raw material and coating . . . . . . . . . . . . . . . . . xi

Figure 1 :

 1 Figure 1: A six-injector vertical-axis Pelton runner [HEN 92].

Figure 2 :

 2 Figure 2: Erosion damage visible on a Pelton turbine bucket at the end of its functional lifetime. ANDRITZ Hydro.

Figure 1 . 1 :

 11 Figure 1.1: Shock front and highly compressed volume on early time of droplet impact [HEY 69].

Figure 1 . 2 :

 12 Figure 1.2: Geometrical construction of wave front with Huygens-Fresnel principle [FIE 99].

Figure 1 . 3 :Figure 1 . 4 :

 1314 Figure 1.3: Time evolution of density during the droplet impact for a perfect fluid, R = 0.1 mm V = 500 m.s -1 and ρ 0 f = 1000 kg.m -3 . Only the x ≥ 0 part is represented due to the symmetric behaviour of the droplet. Corresponding times are (a) 7.98 ns, (b) 18.04 ns, (c) 38.04 ns, (d) 56.96 ns, (e) 95.58 ns, (f) 109.33 ns, (g) 123.24 ns, (h) 138.94 ns, (i) 162.58 ns. [HAL 02a].

Figure 1 . 5 :

 15 Figure 1.5: Water-hammer pressure vs. impact velocity for water droplet on rigid target, with c 0 f = 1647 m.s -1 and ρ 0 f = 1000 kg.m -3 .

Figure 1 . 7 :

 17 Figure 1.7: Spatial distribution of contact pressure for several times after impact (R = 0.1 mm, V = 500 m.s -1 ) [HAL 02b].

1.Figure 1 . 8 :

 18 Figure 1.8: Map of high pressure duration T hp as a function of droplet radius R and impact velocity V for water droplet on rigid target with c 0 f = 1647 m.s -1 .

Figure 1 . 9 :

 19 Figure 1.9: Birth of lateral microjets [HAL 02a].

Figure 1 . 10 :

 110 Figure 1.10: Radial velocity and density distribution on the contact zone when jetting [HAL 02b].

Figure 1 . 11 :

 111 Figure 1.11: Steps of near-side cavitation formation. (a) Born of depression, generated by intercepting shock wave fronts. (b) Cavitation bubbles near to the surface, which can induce erosion when collapsing. [OBR 11].

Figure 1 . 12 :

 112 Figure 1.12: Impact of 10 mm diameter water drop by a metal slider at a velocity of 110 m.s -1 . The onset of jetting is visible in frame c and labelled 'J' is framed d. The shock, labelled 'R', reflects in frame i and focuses, labelled 'F', in frame j. Interframe time is 1 ms. [FIE 89].

Figure 1 . 15 :

 115 Figure 1.15: Example of a large, high-speed, rotating arm-and-spray distributed impact apparatus (courtesy of Bell Aerospace TEXTRON, Buffalo, NY) [AST 10].

1.Figure 1 . 16 :

 116 Figure 1.16: Typical erosion curve [BAK 66].

Figure 1 . 18 :

 118 Figure 1.18: Erosion behaviour of X20Cr13 at different droplet sizes tested at the impact velocity of 488 m.s -1 [AHM 13].

Figure 1 . 19 :

 119 Figure 1.19: Erosion behaviour of 12%Cr steel at different impact velocities for a droplet size φ = 0.2 mm [GER 95].

4.

  These defaults are amplified by fatigue. Finally, erosion is driven by plastic deformation, hardening, intergranular cracks propagation and fatigue mechanism. The cyclic nature of the damage produces a digging by steps (Figure 1.21).

Figure 1 . 20 :

 120 Figure 1.20: Grain ruptures on the initial surface of H300 after 50'000 impacts at impact velocity V = 225 m.s -1 [LUI 13].

Figure 1 . 21 :

 121 Figure 1.21: Fatigue defects observed on H300 after 5 million (on left) and 10 million (on right) impacts, at impact velocity V = 225 m.s -1 [LUI 13].

Figure 1 . 22 :

 122 Figure 1.22: Cross section of eroded stainless steel at different magnification, at impact velocity V = 225 m.s -1 [LUI 13].

Figure 1 . 23 :

 123 Figure 1.23: Intergranular damage observed on Ti-6Al-4V, at impact velocity V = 350 m.s -1 , (a) large scale and (b) high magnification. [KAM 15].

Figure 1 . 24 :

 124 Figure 1.24: Born of triple junctions an microvoids on Ti-6Al-4V, at impact velocity V = 350 m.s -1 , (a) large scale and (b) high magnification. [KAM 15].

  [ROC 72]. Propagation induced by repetition of impacts leads to macro and microcracks [MEN 98], and combination of normal and shear forces (P and S-waves) generates slip and shear bands (Figure 1.26).

Figure 1 . 25 :

 125 Figure 1.25: Grain ejection, striation inside crater and slip bands on Ti-6Al-4V, at impact velocity V = 350 m.s -1 , (a) large scale and (b) high magnification. [KAM 15].

Figure 1 .

 1 Figure 1.26: Stress wave propagation observed on γ-TiAl surface on two regions.[KON 10].

Figure 1 . 27 :

 127 Figure 1.27: Intergranular and intragranular cracking in γ-TiAl, (a) born of cracks, (b) opening of intergranular cracks [KON 10].

Figure 1 . 29 :

 129 Figure 1.29: The two modes of intragranular fractures in γ-TiAl : (a) translamellar crack, (b) interlamellar crack. [KON 10].

Figure 1 . 30 :

 130 Figure 1.30: Areas sensitive to cracking in Pelton buckets [DéL 14].

  .34). The resulting coating has low porosity and high bond strength [KUR 08], and mechanical properties of HVOF coating are checked in [BRA 95]. Figure 1.35 shows the difference of surface state between a non-coated Pelton bucket and two coated ones, after many hours of operation.

Figure 1 .Figure 1 Figure 1 Figure 1 . 35 :

 111135 Figure 1.31: HVOF coating on Pelton buckets. ANDRITZ Hydro.

Figure 2 . 1 :

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

Figure 2 . 2 :

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

Figure 2 . 3 :

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

Figure 2 . 4 :

 24 Figure 2.4: Representation and dimensions of the numerical model.

Figure 2 . 6 :

 26 Figure 2.6: Distribution of probes in 2-D model.

Figure 2 . 7 :

 27 Figure 2.7: Distribution of probes in 3-D model.

Figure 2 . 8 :

 28 Figure 2.8: Impulse distribution of 3-D computation at t = 100 ns after impact for R = 0.5 mm, and V = 100 m.s -1 .

Figure 2 . 9 :

 29 Figure 2.9: Comparison of impulse distribution between 2-D and 3-D models at t = 100 ns after impact for R = 0.5 mm, and V = 100 m.s -1 .

Figure 2 . 10 :Figure 2 . 11 :

 210211 Figure 2.10: Comparison of impulse distribution between full and symmetric models at t = 100 ns after impact for R = 0.5 mm, and V = 100 m.s -1 .

Figure 2 . 12 :

 212 Figure 2.12: Comparison of impulse distribution between full and truncated models at t = 100 ns after impact for R = 0.5 mm, and V = 100 m.s -1 .

Figure 2 . 13 :

 213 Figure 2.13: Comparison of impulse distribution for different particle sizes at t = 64 ns after impact for R = 0.5 mm, and V = 100 m.s -1 .

Figure 2 . 14 :

 214 Figure 2.14: Contact pressure on rigid wall for R = 0.5 mm, and V = 100 m.s -1 .

Figure 2 . 15 :

 215 Figure 2.15: Pressure field p f inside the droplet at t = 150 ns after impact, with a particle size of ∆x = 6 • 10 -7 m.

Figure 2 . 16 :

 216 Figure 2.16: Hydrostatic stress σ H inside the solid at t = 600 ns after impact. The scale compared to the droplet is indicated by the mark x = R. The elements have a size of 5 • 10 -5 m by simplification.

Figure 2 .

 2 Figure 2.17: Pressure inside fluid, hydrostatic and Von Mises stresses (respectively left and right) inside solid for several times after impact : (a) t = 60 ns, (b) t = 260 ns, (c) t = 440 ns. Fluid particles have a size of ∆x = 6 • 10 -6 m. Same isovalue scale for every view, different view scales.

Figure 2 . 18 :

 218 Figure 2.18: Time evolution of hydrostatic and Von Mises stresses at x = 0.1425 mm, z = -2.5 µm. Droplet impacts at t = 130 ns.

Figure 2 . 19 :

 219 Figure 2.19: Time evolution of stress triaxiality ζ at x = 0.1425 mm, z = -2.5 µm.Droplet impacts at t = 130 ns.

Figure 2 .

 2 Figure 2.20: Signed Von Mises Stress σ sign VM inside solid and pressure inside fluid at jetting time T jet ≈ 70 ns after impact. Fluid particles have a size of ∆x = 6 • 10 -6 m.

Figure 2 . 21 : 2 .

 2212 Figure 2.21: Time evolution of principal stresses at x = 0.1425 mm, z = -2.5 µm. Droplet impacts at t = 130 ns.

Figure 2 . 22 :

 222 Figure 2.22: Time evolution of principal stresses ratios r i, j at x = 0.1425 mm, z = -2.5 µm. Droplet impacts at t = 130 ns.

Figure 2 . 23 :

 223 Figure 2.23: Time evolution of principal angle α p at x = 0.1425 mm, z = -2.5 µm.Droplet impacts at t = 130 ns.

Contents 3 . 1

 31 The fatigue phenomenon . . . . . . . . . . . . . . . . . . . . . . . . . . 91 3.1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 3.1.2 Generalities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 3.1.3 Definition of terms used in fatigue . . . . . . . . . . . . . . . . . . 91 3.1.4 The different fatigue domains . . . . . . . . . . . . . . . . . . . . 93 3.1.5 The different fatigue loads . . . .

Figure 3

 3

Figure 3 . 1 :

 31 Figure 3.1: Example of sinusoidal load cycles with R σ < 0.

Figure 3 Figure 3 . 2 :

 332 Figure 3.2: S-N diagram and different fatigue domains.

Figure 3 . 4 :

 34 Figure 3.4: Definition of terms related to normal stress σ h (t).

Figure 3 .

 3 5 shows the decomposition of the shear stress and determination of stress amplitude. The center of circle C is noted O, and the origin of the basis, where (τ uh , τ vh ) = (0, 0) is placed on the material point P. The shear stress τ h (t) is represented by vector -→ PA in Figure 3.5, and is then the combination of alternate stress τ a h (t) (vector -→ OA) and mean stress τ m h (vector -→ PO) (cf. equation (3.15)). Stress amplitude τ am h being the maximum value of the norm of alternate stress τ a h (t) (cf. equation (3.16)), it is equal to the radius of C . The shear range ∆τ h = 2τ am h is sometimes used instead, and is the diameter of C . Some authors define a priori the critical plane, on which the fatigue function is calculated. In those cases, the fatigue function has a more complex expression than the one determining the critical plane. Inside this set of criteria, we can find : Stulen & Cummings [STU 54], Findley [FIN 56], Mc Diarmid [MCD 73, MCD 74], Matake [MAT 77, MAT 80], Mubday & Mitchell [MUN 89] and Galtier & Séguret [GAL 90, GAL 93].

  [YOK 66], Flavenot & Skalli [FLA 82], Dang Van [DAN 71, DAN 84, DAN 89], Froustey [FRO 88], Deperrois [DEP 91] and Robert [ROB 92].

  3.4) : Sines [SIN 55, SIN 59, SIN 81], Marin [MAR 56], Crossland [CRO 56, CRO 70], Dietman & Issler [DIE 74], Kinasoshvili [PIS 76, ALT 96], Kakuno & Kawada [KAK 79], Hashin [HAS 81, YUA 86] Papadopoulos 1 [PAP 87], Altenbach

Figure 3 . 5 :

 35 Figure 3.5: Definition of terms related to tangential stress τ h (t).

Figure 3 . 6 :

 36 Figure 3.6: Endurance domain and two typical load paths G 1 , G 2 . [DAN 89]

Initialization 2 -Figure 3 . 7 :Figure 3

 2373 Figure 3.7: Procedure for erosion simulation.

FatigueFigure 3 . 9 :

 39 Figure 3.9: Procedure for calculation of stress amplitude for the ASME criterion.

Figure 3 . 10 :

 310 Figure 3.10: Procedure for calculation of stress amplitude for the Dang Van 2 criterion.

ConclusionFigure 3 . 11 :

 311 Figure 3.11: Example of mesh with eroded elements and modification of FSI interface (a) Before erosion, (b) After erosion.

Contents 4 . 1

 41 Features and assumptions . . . . . . . . . . . . . . . . . . . . . . . . . . 112 4.2 Results for non-coated surface . . . . . . . . . . . . . . . . . . . . . . . 114 4.2.1 Results of parametric analysis . . . . . . . . . . . . . . . . . . . . 114 4.2.2 Solid mesh convergence checking . . . . . . . . . . . . . . . . . . 120 4.3 Adding the coating layer . . . . . . . . . . . . . . . . . . . . . . . . . . 121 4.3.1 Effects on transient evolution . . . . . . . . . . . . . . . . . . . . . 122 4.3.2 Effects on fatigue resistance and lifetime . . . . . . . . . . . . . . 124 4.3.3 Discussion about results and limitations . . . . . . . . . . . . . . . 127 4.4 Conclusion on prediction of damage due to droplets impacts . . . . . . 127

  Figure 4.2.a and fatigue S-N diagram in Figure 4.2.b.

  .2.a. Its S-N diagram is shown in Figure 4.2.b.

Figure 4 . 1 :

 41 Figure 4.1: Representation and dimensions of numerical model involving coating layer.

Figure 4 . 2 :

 42 Figure 4.2: Raw and coating materials behaviour curves : (a) Traction curve, (b) S-N diagram [BRA 95, WEI 16].

Figure 4 . 3 :

 43 Figure 4.3: Stress amplitude inside solid according to ASME criterion σ am AS on left and Dang Van 2 criterion σ am DV 2 on right, and pressure inside fluid p f at t = 460 ns after impact for the case (φ = 1 mm,V = 100 m.s -1 ). Marks indicate the theoretical position of microjets formation at x = ±R jet . Fluid particles have a size of ∆x = 6 • 10 -6 m.

Figure 4 . 4 :

 44 Figure 4.4: Number of cycles to failure according to ASME criterion N AS on left and Dang Van 2 criterion N DV 2 on right at t = 460 ns after impact for the case (φ = 1 mm,V = 100 m.s -1 ). Marks indicate the theoretical position of microjets formation at x = ±R jet . Fluid particles have a size of ∆x = 6 • 10 -6 m.

Figure 4

 4 Figure4.5: (a) Influence of diameter φ for different velocities V j on the maximum stress amplitude over the mesh according to both criteria, (b) Influence of velocity V for different diameters φ i on the maximum stress amplitude over the mesh according to both criteria.

Figure 4

 4 Figure 4.6: (a) Influence of diameter φ for different velocities V j on minimum number of cycles to failure over the mesh according to both criteria, (b) Influence of velocity V for different diameters φ i on minimum number of cycles to failure over the mesh according to both criteria.

Figure 4

 4 Figure 4.7: (a) Influence of diameter φ for different velocities V j on number of eroded elements according to both criteria, (b) Influence of velocity V for different diameters φ i on number of eroded elements according to both criteria. Both graphics are for N lim = 5 • 10 7 cycles.

Figure 4

 4 Figure 4.8: (a) Influence of diameter φ for different velocities V j on erosion depth according to both criteria, (b) Influence of velocity V for different diameters φ i on erosion depth according to both criteria. Both graphics are for N lim = 5 • 10 7 cycles.

Figure 4 . 9 :

 49 Figure 4.9: Fatigue function and non-eroded elements for N lim = 5 • 10 7 cycles according to ASME criterion E AS on left and Dang Van 2 criterion E DV 2 on right at t = 460 ns after impact for the case (φ = 0.5 mm,V = 100 m.s -1 ). Marks indicate the theoretical position of microjets formation at x = ±R jet . Fluid particles have a size of ∆x = 6 • 10 -6 m.

Figure 4 . 10 :

 410 Figure 4.10: Fatigue function and non-eroded elements for N lim = 5•10 7 cycles according to ASME criterion E AS on left and Dang Van 2 criterion E DV 2 on right at t = 460 ns after impact for the case (φ = 0.5 mm,V = 200 m.s -1 ). Marks indicate the theoretical position of microjets formation at x = ±R jet . Fluid particles have a size of ∆x = 6 • 10 -6 m.

Figure 4 . 11 :

 411 Figure 4.11: Fatigue function and non-eroded elements for N lim = 5•10 7 cycles according to ASME criterion E AS on left and Dang Van 2 criterion E DV 2 on right at t = 460 ns after impact for the case (φ = 1 mm,V = 100 m.s -1 ). Marks indicate the theoretical position of microjets formation at x = ±R jet . Fluid particles have a size of ∆x = 6 • 10 -6 m.

  and 4.18 respectively, and the comparison is shown for both models : non-coated on the left and coated on the right. Isovalues scales are reduced in order to highlight the stresses gradients.The material data of coating given in section 4.1 and equation (1.1) from Chapter 1 allow calculating the velocities of P and S-waves inside the coating : c coat L = 7238 m.s -1 and c coat T

  and 4.18 represent typical times for the observation of wave propagation :

Figure 4 . 12 :

 412 Figure 4.12: Fatigue function and non-eroded elements for N lim = 5•10 7 cycles according to ASME criterion E AS on left and Dang Van 2 criterion E DV 2 on right at t = 460 ns after impact for the case (φ = 1 mm,V = 200 m.s -1 ). Marks indicate the theoretical position of microjets formation at x = ±R jet . Fluid particles have a size of ∆x = 6 • 10 -6 m.

Figure 4 . 13 :

 413 Figure 4.13: Fatigue function and non-eroded elements for N lim = 5•10 7 cycles according to ASME criterion E AS on left and Dang Van 2 criterion E DV 2 on right at t = 460 ns after impact for the case (φ = 2 mm,V = 100 m.s -1 ). Marks indicate the theoretical position of microjets formation at x = ±R jet . Fluid particles have a size of ∆x = 6 • 10 -6 m.

Figure 4 . 14 :

 414 Figure 4.14: Fatigue function and non-eroded elements for N lim = 5•10 7 cycles according to ASME criterion E AS on left and Dang Van 2 criterion E DV 2 on right at t = 460 ns after impact for the case (φ = 2 mm,V = 200 m.s -1 ). Marks indicate the theoretical position of microjets formation at x = ±R jet . Fluid particles have a size of ∆x = 6 • 10 -6 m.

Figure 4 .

 4 Figure 4.15: Dang Van 2 stress amplitude σ am DV 2 for standard mesh (left) and refined one (right) in the highest stress region, at t = 460 ns after impact for (φ = 1 mm,V = 100 m.s -1 ).Fluid particles have a size of ∆x = 6 • 10 -6 m.

Figure 4 . 17 :

 417 Figure 4.17: Comparison of hydrostatic stress σ H between the pure basic material (left) and the coated one (right) at different frames for the case (φ = 1 mm,V = 100 m.s -1 ) : (a) t = 50 ns, (b) t = 90 ns, (c) t = 170 ns, (d) t = 360 ns after impact. Fluid particles have a size of ∆x = 6 • 10 -6 m.

Figure 4 . 18 :Figure 4 . 19 :Figure 4 . 20 :

 418419420 Figure 4.18: Comparison of Von Mises stress σ VM between the pure basic material (left) and the coated one (right) at different frames for the case (φ = 1 mm,V = 100 m.s -1 ) : (a) t = 50 ns, (b) t = 90 ns, (c) t = 170 ns, (d) t = 360 ns after impact. Fluid particles have a size of ∆x = 6 • 10 -6 m.

Figure 4 . 22 :

 422 Figure 4.22: Comparison of Dang Van 2 stress amplitude σ am DV 2 between pure basic material (left) and coated one (right) at t = 460 ns after impact for φ = 1 mm and V = 100 m.s -1 . Marks indicate the theoretical position of microjets formation at x = ±R jet . Fluid particles have a size of ∆x = 6 • 10 -6 m.

Figure 4 . 23 :

 423 Figure 4.23: Comparison of number of cycles to failure according to ASME criterion N AS between pure basic material (left) and coated one (right) at t = 460 ns after impact for φ = 1 mm and V = 100 m.s -1 . Marks indicate the theoretical position of microjets formation at x = ±R jet . Fluid particles have a size of ∆x = 6 • 10 -6 m.

Figure 4 . 24 :

 424 Figure 4.24: Comparison of number of cycles to failure according to Dang Van 2 criterion N DV 2 between pure basic material (left) and coated one (right) at t = 460 ns after impact for φ = 1 mm and V = 100 m.s -1 . Marks indicate the theoretical position of microjets formation at x = ±R jet . Fluid particles have a size of ∆x = 6 • 10 -6 m.

Figure 4 . 25 :

 425 Figure 4.25: Comparison of fatigue function for N lim = 3 • 10 9 cycles according to ASME criterion E AS between pure basic material (left) and coated one (right) at t = 460 ns after impact for φ = 1 mm and V = 100 m.s -1 . Marks indicate the theoretical position of microjets formation at x = ±R jet . Fluid particles have a size of ∆x = 6 • 10 -6 m.

Figure 4 . 26 :

 426 Figure 4.26: Comparison of fatigue function for N lim = 3 • 10 9 cycles according to Dang Van 2 criterion E DV 2 between pure basic material (left) and coated one (right) at t = 460 ns after impact for φ = 1 mm and V = 100 m.s -1 . Marks indicate the theoretical position of microjets formation at x = ±R jet . Fluid particles have a size of ∆x = 6 • 10 -6 m.

d

  Compression waves velocity inside fluid when jetting m.s -1 c L Elastic compression waves velocity of solid m.s -1 c coat L Elastic compression waves velocity of coating m.s -1 c T Elastic shear waves velocity of solid m.s -1 c coat T Elastic shear waves velocity of coating m.s -1 Dimension

Y

  related to the plane normal to h Pa σ am h Normal stress amplitude related to the plane normal to h Pa σ m h Mean normal stress related to the plane normal to h Pa σ max h Maximum normal stress related to the plane normal to h Pa Initial yield strength of solid Pa σ 0,coat Y Initial yield strength of coating Pa σ h Normal stress related to the plane normal to h Pa ς Surface tension N.m -1 τ -1 Endurance limit under symmetrical alternate torsion Pa τ oct Tangential octahedral stress Pa τ uh Component of τ h on u h Pa τ vh Component of τ h on v h Pa τ a Alternate shear Pa τ h Tangential stress related to the plane normal to h Pa τ a h Alternate tangential stress related to the plane normal to h Pa τ am h Tangential stress amplitude related to the plane normal to h Pa

Figure 27 :

 27 Figure 27: Front d'onde et volume comprimé en début d'impact [HEY 69].

Finalement

  , l'érosion est pilotée par la déformation plastique, l'écrouissage, la propagation de fissures intergranuilaires et la fatigue mécanique. La nature cyclique du dommage est confirmée par une observation de la zone d'érosion qui prend une forme d'escalier, visible sur la Figure29.

Figure 29 :

 29 Figure 29: Dommage par fatigue observé sur de l'acier H300 après 5 millions d'impacts (gauche) et 10 millions (droite) à vitesse V = 225 m.s -1 [LUI 13].

Figure 30 :

 30 Figure 30: Représentation et dimensions du modèle numérique.

Figure 31 :

 31 Figure 31: Pression de contact sur le mur rigide pour R = 0.5 mm et V = 100 m.s -1 .

Figure 32 :Figure 33 :

 3233 Figure 32: Pression dans le fluide, contraintes hydrostatique et de Mon Mises (respectivement à gauche et à droite) dans le solide en différents instants après impact : (a) t = 60 ns, (b) t = 260 ns, (c) t = 440 ns. Les échelles d'isovaleurs sont identiques, mais pas la vue.

InitialisationFigure 34 :

 34 Figure 34: Procédure de la simulation d'érosion.

2

 2 , l'étendue de contrainte ∆σ est obtenue comme étant la plus grande contrainte de Rankine de la différence de contraintes ∆σ (m,n) , et l'amplitude de contrainte en est la moitié : Critère de DangVan 2 [DAN 89]

Figure 35 :

 35 Figure 35: Fonction de fatigue et éléments non-érodés pour N lim = 5 • 10 7 cycles d'après le critère ASME E AS (gauche) et celui de Dang Van 2 E DV 2 (droite), au temps t = 460 ns après impact, pour le cas (φ = 1 mm,V = 100 m.s -1 ). Les cotes indiquent la position théorique de formation des microjets en x = ±R jet .

  Tout d'abord, le mécanisme d'érosion a été décrit en deux parties. Les phénomènes de propagation d'ondes dans la goutte permettent d'expliquer la formation des microjets et

Figure 36 :

 36 Figure 36: Comparaison de la contrainte de Von Mises σ VM entre le modèle en acier pur (gauche) et celui avec coating (droite) en différents instants pour le cas (φ = 1 mm,V = 100 m.s -1 ) : (a) t = 50 ns, (b) t = 90 ns, (c) t = 170 ns, (d) t = 360 ns après impact.

Figure 37 :

 37 Figure 37: Comparaison de la fonction de fatigue pour N lim = 3 • 10 9 cycles d'après le critère de Dang Van 2 E DV 2 entre le modèle en acier pur (gauche) et celui avec coating (droite) à t = 460 ns après impact pour φ = 1 mm et V = 100 m.s -1 .

  FOLIO ADMINISTRATIF THESE DE L'UNIVERSITE DE LYON OPEREE AU SEIN DE L'INSA LYON NOM : COUDOUEL DATE de SOUTENANCE : 26/10/17 Prénoms : Guillaume, Jean, André TITRE : Toward a numerical predictive method based on fatigue analysis for droplet impingement erosion NATURE : Doctorat Numéro d'ordre : 2017LYSEI101 Ecole doctorale : MEGA Spécialité : Génie Mécanique RESUME : Le but du travail présenté est la compréhension puis la simulation numérique des mécanismes d'érosion des augets de turbine Pelton par impacts répétés de gouttes d'eau dans le but de prédire la durée de vie des composants. Tout d'abord, les phénomènes de propagation d'ondes dans les milieux fluide et solide sont étudiés. Cela permet de mettre en lumière l'évolution temporelle et la distribution spatiale des pressions de contact, et l'apparition de microjets par éjection supersonique du fluide au contact. Les études expérimentales de l'érosion par gouttes d'eau traduisent un dommage basé sur la fissuration par fatigue. Des simulations numériques en dynamique rapide couplées fluide-structure sont alors effectuées. Le domaine solide est discrétisé par la Méthode des Éléments Finis (MEF), et le domaine fluide par la méthode « Smoothed Particle Hydrodynamics » (SPH), qui est une méthode particulaire (sans maillage) particulièrement adaptée aux grandes distorsions et au suivi des surfaces libres. L'analyse des états de contraintes vient corroborer la nature cyclique de l'endommagement. La simulation d'érosion est alors réalisée à l'aide de critères de fatigue multiaxiaux. Le choix se porte vers un premier critère général de l'American Society of Mechanical Engineers (ASME), utilisant les valeurs principales des différences de contraintes au cours du temps. Le second choix concerne un critère à plan critique : le critère de Dang Van 2. Il traite séparément la contrainte hydrostatique et le cisaillement alterné maximal local. Ces critères permettent de définir les régions érodées du solide au bout d'un nombre d'impact donné, ce qui fait de cette démarche une méthode prédictive. Une étude paramétrique pour différentes tailles de gouttes et vitesses d'impact est ensuite réalisée, puis on évalue l'influence de la présence d'une couche de coating. MOTS-CLÉS : Impact de goutte, Érosion, Couplage fluide-structure, Fatigue, Éléments Finis, Smoothed Particle Hydrodynamics. Laboratoire de recherche : Laboratoire de Mécanique des Contacts et des Structures (LaMCoS) Directeur de thèse : Alain COMBESCURE Présidente de jury : Catherine GARDIN Composition du jury : Alain COMBESCURE, Vincent FAUCHER, Catherine GARDIN, Habibou MAITOURNAM, Jean-Christophe MARONGIU, Emmanuelle SALLÉ

  Conclusion on the mechanism of droplet impact erosion . . . . . . . . . . 2 Numerical methods for droplet impact 2.1 Strong form of continuum mechanics . . . . . . . . . . . . . . . . . . . . 2.1.1 Different configurations in continuum mechanics . . . . . . . . . 2.1.2 Conservations laws in continuum mechanics . . . . . . . . . . . 2.1.3 The Solid sub-domain . . . . . . . . . . . . . . . . . . . . . . . 2.1.4 The fluid sub-domain . . . . . . . . . . . . . . . . . . . . . . . . 2.1.5 Coupling conditions . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Weak form and numerical discretization . . . . . . . . . . . . . . . . . . 2.2.1 FEM for solid sub-domain . . . . . . . . . . . . . . . . . . . . . The fatigue phenomenon . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.2 Generalities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.3 Definition of terms used in fatigue . . . . . . . . . . . . . . . . . 3.1.4 The different fatigue domains . . . . . . . . . . . . . . . . . . . 3.1.5 The different fatigue loads . . . . . . . . . . . . . . . . . . . . . 3.1.6 Conclusion on fatigue phenomenon . . . . . . . . . . . . . . . . 3.2 Multiaxial fatigue criteria . . . . . . . . . . . . . . . . . . . . . . . . . . Conclusion on fatigue developments . . . . . . . . . . . . . . . . . . . .

	3 Fatigue developments
	3.1
	i

1.2.1 Test rigs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.2 Macroscopic scale . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.3 Material science based models . . . . . . . . . . . . . . . . . . . 1.2.4 Coating . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 2.2.2 SPH for fluid sub-domain . . . . . . . . . . . . . . . . . . . . . 2.2.3 Time integration . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 General ASME criterion . . . . . . . . . . . . . . . . . . . . . . 3.2.2 Specific criteria . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Choice of fatigue criterion . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.1 Presentation of Dang Van criterion [DAN 71, DAN 84] . . . . . . 3.3.2 Determination of the critical plane . . . . . . . . . . . . . . . . . 3.3.3 Second version of Dang Van criterion [DAN 89] . . . . . . . . . 3.4 Fatigue erosion predicting tool . . . . . . . . . . . . . . . . . . . . . . . 3.4.1 General procedure . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.2 Protocol for stress amplitude calculation . . . . . . . . . . . . . . 3.4.3 Procedure for new FSI boundary determination . . . . . . . . . . 3.5

  After a larger number of impacts, neighbour grains support the same damage mechanism and are ejected or fractured (step 2.). Microckacks are intergranular type, which impairs the surface condition, and move in parallel to the surface and propagate in depth. The damage zone can be larger than the droplet itself (Figures 1.21 and 1.22).

	1. Physical phenomenon of droplet impact erosion
	water-hammer pressure and the additionnal loads from cavity collapse [FIE 12].
	1.2.3 Material science based models
	A first mesoscopic description of erosion mechanism is proposed by Luiset [LUI 13] as
	the following steps :
	1. The first impacts start to erode grain boundaries and generate pits between grains.
	Then, microcracks appear at the bottom of these pits.
	2. Next, material is removed from surface by two damage modes : a) grain ejection
	which can produce triple junctions, b) grain fracture (Figure 1.20).
	3.
	1.2.2.5 Discussion on phenomenological empirical models
	All these models are phenomenologycal ones in the sense that they are based on fitting
	of experimental results. They do not consider the liquid properties, such as density or
	viscosity. Erosion rate varies at 2 nd to 2.5 th power with liquid density and 1.3 th to 2 nd
	power with the liquid viscosity. An increase in temperature of the impacting liquid gene-
	rally increases the erosion slightly. This effect results from the increased shear damage of
	the surface caused by the evolving lateral jet flow [BOW 61, HAN 66, HEY 92, SMI 66].
	It has been argued that this interaction between the lateral outflow jetting and the sur-
	face discontinuities makes a significant contribution to the damage of lightly roughened
	surfaces [ADL 79, HUA 12]. According to Ma and Mann [MA 15, MAN 02], effects of
	microjets are significant after first material removals. An investigation of erosion acting
	on low-pressure side of blades in steam turbines was produced by Azevedo and Staniša
	[AZE 09, STA 95] and concluded that erosion comes from fatigue cracking. We can note
	that damage threshold value should not be affected by water droplet impact for brittle
	materials, because they are strong in compression and should not fail when loaded by

  × 1] gradient of shape function B i from equation (2.69), becomes the [6 × 3] matrix B i . This allows writing the virtual strain rate and strain (2.77) as :

				2.78)
	with
	{σ} =	σ 11 σ 22 σ 33 σ 12 σ 23 σ 31	(2.79)
	and
	{ε} =	ε 11 ε 22 ε 33 2ε 12 2ε 23 2ε 31	(2.80)
	Similarly, the virtual strain tensor δD from equation (2.59) can also be expressed with
	Voigt notation method :
	{δD} =	δD 11 δD 22 δD 33 2δD 12 2δD 23 2δD 31	(2.81)
	Then, the [3

  2. Numerical methods for droplet impact (2.98) (like the shape functions in Finite Element Method). But this expression works for the fluid particles whose support domain is not truncated by the solid boundary or free surface. More informations about the truncature problem are given by Li and Marongiu [LI 13, MAR 07].

	2.2.2.2 SPH discretization of fluid sub-domain
	The most common SPH method in fluid mechanics uses Lagrangian description and is
	denoted by Classic Monhagan SPH formulation [MON 92, MON 05]. But in this work,
	the equations driving the behaviour of the fluid and the gradients are preferred in their
	Arbitrary Eulerian-Lagrangian (ALE) form (section 2.1.1) for conservation purposes (cf.
	equation 2.43) :

Table 2 . 1 :

 21 Material data for numerical simulations.

	Solid	Fluid

Table 2 . 2 :

 22 CPU time for different 2-D fluid models.

			Model type (2-D)
	Symmetry	No	Yes	Yes	Yes	Yes
	Truncation height	-	-	R	R/2	R/2
	Particle size [m]	6 • 10 -7 6 • 10 -7 6 • 10 -7 6 • 10 -7 3 • 10 -7
	CPU time [s]	279 000 75 500 48 600 18 600 105 700

  Fatigue is the modification of material properties under changing and repeated loads or deformations. This term describes mainly the failures that leads to cracking and material fracture. Most of critical mechanical components break under alternate loads, like crankshaft, connecting rods, pinions in automotive industry, landing gear, turbines blades and compressor in aerospace industry, or turbine buckets in energy production. Among damage modes, 50% to 90% of failures are caused by fatigue [LU 02]. Scientists have studied the fatigue damage since the XIX th century during the Industrial Revolution, and the specificity of fatigue is that loads can be low compared to fracture resistance or yield strength. Fatigue damage can be caused by several types of loads, the most basic being an external alternate load. If residual stresses are located in the material, they are superposed to external loads. These loads can be combined with other kinds of loads, like creep or thermal fatigue, caused by temperature variation. A corrosive environment accelerates the failure of the material.

	. . . . . . . . . . . . . . . . . . 94
	3.1.6 Conclusion on fatigue phenomenon . . . . . . . . . . . . . . . . . 95
	3.2 Multiaxial fatigue criteria . . . . . . . . . . . . . . . . . . . . . . . . . . 95
	3.2.1 General ASME criterion . . . . . . . . . . . . . . . . . . . . . . . 95

3.2.2 Specific criteria . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 3.3 Choice of fatigue criterion . . . . . . . . . . . . . . . . . . . . . . . . . . 100 3.1.1 Introduction

  The oldest are from Hohenemser & Prager [HOH 33] and Gough & Pollard [GOU 35, GOU 52]. The most recent criteria are proposed by Lee [LEE 81, LEE 85, LEE 89, YOU 96].

Table 4 . 1 :

 41 Material data for raw material and coating.

		Raw material		Coating
	ρ 0 s	7700 kg.m -3 ρ 0,coat s	13500 kg.m -3
	E 0 200 GPa	E coat 0	600 GPa
	ν	0.288	ν coat	0.24

  Les turbines hydrauliques peuvent connaître de sévères dommages durant leur fonctionnement, ceci à cause de la mauvaise qualité de l'eau ou de ses conditions d'écoulement. Ces déteriorations sont responsables de perte de puissance et de l'augmentation des coûts de maintenance, ainsi que de la mise en danger de la sécurité des installations. Il est donc intéréssant pour les fabricants et exploitants de centrales hydrauliques de réduire l'intensité du dommage et protéger les composants par traitement de surface pour améliorer leur durée de vie. Cependant, des outils de prédiction fiables et précis de ce dommage sont à ce jour inexistants. Le travail présenté ici concerne l'érosion causée par les impacts répétés de gouttes d'eau sur les augets de turbine Pelton. En effet, pour les turbines Pelton à haute chute (> 1000 m), les jets d'eau impactant les augets sont composés d'un corps principal liquide entouré de gouttelettes. Les observations actuelles montrent que les régions exposées aux impacts de gouttes sont sujettes à un type d'érosion particulier. Premièrement, le mécanisme d'érosion est présenté, puis la simulation numérique du dommage est réalisée par le biais d'un couplage fluide-structure et de critères de fatigue.Cette partie présente le mécanisme responsable de l'érosion des structures métalliques par impact de gouttes. Tout d'abord, les phénomènes se produisant dans la goutte sont décrits, comme la pression "coup de bélier" et la propagation d'ondes. La pression maximale au contact apparait lorsque les jets latéraux emergent de la goutte sur la surface de contact. La chute de pression dans la zone opposée au point de contact peut généner un phenomène de cavitation. Ensuite, les mécanismes d'érosion en eux-mêmes dans le solide sont présentés. Ils sont répartis en deux groupes distincts : les mécanismes décrits à l'échelle macroscopique qui sont issus de modèles phénoménologiques d'une part, et les mécanismes à l'échelle mésoscopique décrivant la propagation de fissures intergranulaires aboutissant à la rotation puis à l'éjection des grains, ce qui produit des cupules en surface, obervables macroscopiquement. HEY 69,LI 11], si une goutte à faible diamètre impacte une structure solide à vitesse élevée, les effets visqueux et de tensition superficielle peuvent être négligés devant ceux de nature inertielle. En effet, une goutte de rayon R = 0.1 mm à vitesse V = 500 m.s -1 procure les nombres de Reynolds de de Weber Re = 50 000 et We = 350 000. De plus, les résultats numériques de [HAL 02b] montrent l'absence d'effets de convection thermique dans le comportement du fluide. Toutes ces simplifications permettent de modéliser la réponse dynamique de la goutte par les équations d'Euler dans les simulations numériques.

	0.1.1 Physique de l'impact de goutte
	Parameter Description Parameter Description 0.1.1.1 Propagation des ondes	Unit (ISO)	Unit (ISO)
	X	s m	Material or Lagrangian coordinates Mean deviatoric stress	Pa	m
	x D'après [HAL 02b, Spatial or Eulerian coordinates S Piola-Kirchhoff 2 stress	Pa	m
	δv τ h	δD δL	Virtual velocity Tangential stress vector related to the plane normal to h Virtual strain rate Virtual velocity gradient Résumé étendu m.s -1 Pa s -1 s -1
	τ a h	δW	Alternate tangential stress vector related to the plane normal to h Virtual rotating rate s -1	Pa
	τ m h χ 0.0 Introduction Mean tangential stress vector related to the plane normal to h Reference coordinates Strain field inside solid -ε s ε e s Elastic part of strain field -	Pa m
	{ε}	ε	p s	Infinitesimal strain in Voigt notation Plastic part of strain field	-	-
	{σ}			εp	Cauchy stress in Voigt notation Plastic strain rate	s -1	Pa
				σ		Cauchy stress tensor	Pa
		σ f		Cauchy stress field inside fluid	Pa
		σ s		Cauchy stress field inside solid	Pa
	Matrices σ (m)		Cauchy stress at time step t m	Pa
	∆σ (m,n)	Stress difference between time steps t m and t n	Pa
	Parameter Description	Unit (ISO)
				1		Identity matrix	-
	1 d B i Sets and domains Unit matrix of the dimension of the problem d Gradient of shape functions associated to node i	--
				C		Elastic stiffness tensor (4 th order tensor)	Pa
	[C] Parameter Description Elastic stiffness matrix in Voigt notation	Pa
	B p	e		Deviatoric strain Principal basis	-
	B r		E		Green-Lagrange strain Reference basis	-
	C		F		Jacobian or gradient of motion Smallest circle circumscribed to T	-
	H ep D i		Elastoplastic tangent modulus (4 th order tensor) Influence area of the kernel function centred on particle i	Pa
	ω max ∂D i K ω σ L G Ω crit M T Ω r N Γ I P Γ s	Highest natural frequency of mesh Stiffness matrix Area boundary of the kernel function centred on particle i N.m -1 Signal pulse of stress load Velocity gradient s -1 Load path in (σ H , τ a h ) space Critical frequency for undamped system Mass matrix kg Load path in tangential stresses space Rotational speed of Pelton runner Nominal stress Pa Fluid-structure interface Piola-Kirchhoff 1 stress Pa Solid boundary	s -1 rad.s -1 s -1 rad.s -1
	Γ t		s		Deviatoric stress Load-constrained boundary	Pa
	Γ v	s a		Alternate deviatoric stress Kinematic-constrained boundary	Pa

0.1 Explication physique de l'érosion par choc de goutte

  OBR 11], une onde de choc traversant un milieu liquide est susceptible de causer l'apparition de cavitation. En effet, [FIE 12] affirme que la superposition des ondes de relaxation en sommet de goutte provoque une chute locale de pression et donc des contraintes de tension (traction), ce qui est connu pour générer de la cavitation. Plus la vitesse d'impact V est élevée, plus le risque de cavitation est grand. Cependant, la cavitation n'est pas considérée dans les simulations numériques de ce travail.En travaillant sur l'érosion des aubes de turbines à vapeur, [BAK 66] propose une description temporelle du mécanisme de d'endommagement macroscopique en quatre étapes : (i) La première phase est généralement nommée période d'incubation car il n'y a pas d'enlèvement significatif de matière mais seulement une modification de l'état de surface.

	0.1.1.3 Cavitation
	D'après [0.1.2 Mécanisme d'érosion
	0.1.2.1 Modèle phénoménologique
	front d'onde se progageant (iv) Finalement, le taux d'érosion diminue puis reste constant, et peut s'annuler dans
	certains cas.	surface libre de la goutte
	0.1.2.2 Modèle basé sur la science des matériaux
	Une description mésoscopique du mécanisme d'érosion est proposée par [KAM 15,
	c f LUI 13] et consiste dans les étapes suivantes :
	1. Les premiers impacts afaiblissent les joints de grains et creusent des puits à la place.
	microjet 2. Ensuite, l'enlèvement de matériau commence par deux modes : a) éjection des
	grains qui engendre des points triples, b) rupture des grains (translamellaire ou
	interlamellaire).
	3. Après un plus grand nombre d'impacts, les grains périphériques sont sujets au
	Figure 28: Formation des microjets latéraux [HAL 02a]. même mode d'endommagement et sont soit éjectés, soit fracturés (étape 2.). Les
	fissures sont de type intergranulaire, ce qui dégrade l'état de surface, et se propa-
	La vitesse d'éjection du fluide peut être très supérieure à la vitesse d'impact V , voire f . Cependant, il n'existe pas à ce jour d'expression analytique de à la vitesse initiale du son c 0 gent en profondeur parallèlement à la surface. Des cupules aparaissent en surface,
	la pression maximale, mais [HAL 02b] a proposé une méthode pour déterminer le temps
	d'apparition du jet T jet .

(ii) Durant la deuxième étape, le taux d'érosion augmente de manière constante.

(iii) Le taux d'érosion atteint sa valeur maximale et se stabilise. et la zone de dommage peut être bien plus large que la goutte elle-même. 4. Tous ces défauts sont amplifiés par fatigue.

  Cette section traite de la simulation transitoire de la goutte d'eau impactant le corps solide. Ensuite, un calcul de fatigue dans le solide par post-traitement permet d'estimer l'intensité du dommage et de le quantifier dans le temps, et donc de prédire la durée de vie du composant. Le domaine solide est modélisé par la Méthode des Élements Finis (MEF) grâce au code de dynamique explicite EuroPlexus R [CEA 02], qui est développé conjointement par le Commissariat à l' Énergie Atomique et aux Énergies Alternatives (CEA) et par le Centre Commun de Recherche de la Commission Européenne (EC/JRC). Ce code est adapté à la dynamique rapide fortement non-linéaire avec prise en compte de l'érosion. Pour le domaine fluide, le choix se porte sur le code ASPHODEL développé en interne par ANDRITZ Hydro. Ce code utilise la méthode particulaire (sans maillage) Smoothed Particle Hydrodynamics (SPH), qui est particulièrement adapté au traitement des grandes distorsions, ASPHODEL gérant bien le suivi des surfaces libres. L'intéraction fluidestructure est assurée par le code de couplage fort développé par [NU Ñ 17b, NU Ñ 17a],

0.2 Simulation numérique de l'érosion par choc de goutte qui a l'avantage d'être énergétiquement conservatif à l'interface pour des pas de temps identiques entre domaines fluide et solide.

  2. Les fonctions de fatigue sont alors mises sous la forme :

	E AS (N lim ) =	σ am AS σ lim	E DV 2 (N lim ) =	σ am DV 2 σ lim	(12)
	0.3.2.1 Critère général ASME [MAL 08]

  max t∈[0,T c ] {τ a (t) + α DV σ H (t)} (17) Le cisaillement alterné τ a (t) s'obtient de la manière suivante : À partir du déviateur des contraintes s(t) = σ(t)σ H (t)1, on définit la contrainte déviatorique alternée s a (t) : s a (t) = s(t) -s m (18) où s m repsésente l'état moyen du déviateur des contraintes au cours du cycle de chargement (centre de la plus petite hypersphère circonscrite au trajet de charge dans l'espace déviatorique). Le cisaillement alterné est obtenu par application du critère de Tresca aux contraintes déviatoriques alternées s a (t), ce qui revient à calculer le plus grand des cisaillements maximaux. Ces derniers se situent sur les plans orientés à 45 • des plans contenant les axes principaux des contraintes. De plus on sait que le déviateur s(t) possède les mêmes axes principaux que σ(t), et le critère de Tresca est le même pour les deux contraintes. Le cisaillement alterné est donc : τ a (t) = Tresca {s a (t)} = max III) sont donnés en fonction des valeurs propres de s a (t) dénotées s a K (t) :

						K	{τ a K (t)}	(19)
	où les cisaillements principaux τ a
	τ a I (t) =	|s a II (t) -s a III (t)| 2	τ a II (t) =	|s a III (t) -s a I (t)| 2	τ a III (t) =	|s a I (t) -s a II (t)| 2	(20)
	Finalement le cisaillement alterné peut être exprimé comme suit :
		τ a (t) =	1 2	max {|s a I (t) -s a II (t)| , |s a II (t) -s a III (t)| , |s a III (t) -s a I (t)|}	(21)
	0.3.3 Étude paramétrique du modèle sans coating

K (t) (K = I, ...,

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

For traction test, the negative values are avoided because it leads to compression, and eventually to buckling, which ruins the test specimen. Moreover, the sample behaviour is generally different in traction and compression, stress-strain curve can be non-antisymetric : σ(-ε) = -σ(ε).

As the simulations are performed in

2-D for fluid, and plane strains for solid, the droplet is modelled as a cylinder, thus m L = ρ 0 f πφ 2 /4.

Remerciements