Hybrid Smoothed-Particle Hydrodynamics/Finite Element Method Simulation of Water Droplet Erosion on Ductile Metallic Targets

: Erosion of metallic surfaces due to the permanent impact of high-speed water droplets is a signiﬁcant concern in diverse industrial applications like turbine blades, among others. In the initial stage of water droplet erosion, there is an incubation regime with negligible mass loss whose duration is strongly dependent on water droplet sizes and velocities, the initial state of the surface, and the material properties of the target. The prediction of the incubation period duration is one of the main topics of research in the ﬁeld. In this work, the interaction of the water droplets with a metallic surface is simulated using a hybrid Smoothed-Particle Hydrodynamics/Finite Element Method modeling scheme. The effect of multiple random impacts on representative target areas for certain ranges of impact angles and velocities was studied using a combination of simple material and damage models for the target surface of Ti-6Al-4V titanium alloy. The simulation is able to reproduce the main dependencies of the incubation regime and the ﬁrst stages of water droplet erosion on the impact angle and velocity as reported in the literature. This framework can be considered a foundation for more advanced models with the goal of a better understanding of the physical mechanisms behind the incubation regime in order to devise strategies to extend it in real applications


Introduction
The repeated impact of high-velocity water droplets on metallic surfaces produces an accumulation of surface and bulk damage that eventually evolves into erosion [1][2][3].Water droplet erosion is a serious problem in low-pressure steam turbines [4][5][6], aeroengine fan blades [7], and wind turbine blades [8][9][10].During water droplet impact, high amplitude stress waves propagate inside the target material eventually surpassing the dynamic fracture strength of the target in some locations, in particular those with microstructural discontinuities [11].The repeated effect of those stress waves has also been treated as an analog of fatigue failure [12].Finally, resulting plastic deformation, damage, and microcracks are progressively accumulated eventually leading to the onset of mass loss.
In the initial stage of water droplet erosion (WDE), there is an incubation period with negligible mass loss whose duration is strongly dependent on water droplet size and velocity distributions, the initial state of the target surface, and material properties of the target [13][14][15].Material properties like yield, tensile, and fracture strength affect the ability of materials to resist erosion [3,16].Hardness is directly related to the surface resistance to plastic deformation and definitively is a main factor for determining the erosion resistance of specific material alloys [1,17].The duration of the incubation period is also dependent on the degree of strain hardening occurring after multiple impacts [15,18].Clearly, there is not an intrinsic mechanical property that individually accounts for water droplet erosion and most of the experimental studies consider a combination of a number of mechanical properties [17].From the point of view of industrial applications, the prediction of the incubation period duration and formulation of strategies for its extension is one of the main topics of research [3].
It is important to emphasize that water droplet erosion occurs only above a certain threshold impact velocity for a fixed average droplet size.Namely, for velocities below this threshold, the impact energy is not sufficient for damage accumulation even for a large number of impacts and, therefore, no significant surface deterioration nor material removal will ever occur [19].Ibrahim and Medraj recently proposed an analytical model to predict this threshold velocity using material properties and impact conditions [20].This model is based on the dynamical threshold for crack propagation and other material resistance against erosion parameters that can be measured in rotation erosion test rigs.
The complex interplay between fluid and solid mechanics that drive abrasive wear of surfaces exposed to the impact of water droplets represents a major theoretical and experimental challenge.For example, in order to capture, with a high-speed camera, the dynamical interaction of a single droplet with a metallic surface, the temporal resolution is in the order of the magnitude of billions per second [21].Even after decades of extensive analytical and experimental studies based on empirical modeling, a better understanding of precise physical mechanisms responsible for fracture under erosion conditions is still missing.Diverse numerical simulation methods have been applied to complement the analysis of the impact of droplets on solid surfaces [22,23].The major challenge is to couple the fluid dynamics of the droplets, the impact, and the mechanics of damage and fracture of materials within a coherent numerical simulation framework.
Concerning recent works on numerical simulation, Chidambaram et al. solved numerically the three-dimensional compressible Navier-Stokes equations for the fluid droplets flowing over the blade of a gas turbine engine [24].An empirical model of the resulting erosion rate was coupled to the simulation and results were compared with experiments.Castorini et al. applied an analogous simulation approach for rain erosion in wind turbine blades [25].Armizadeh et al. proposed a computational framework for rain erosion prediction in wind turbine blades by combining probability theory and stochastic analysis for rain droplet size and velocity distributions, computational fluid dynamics, transient finite element modeling, and fatigue analysis [26,27].
Simulating the impact and splashing of water droplets interacting with an already eroded surface requires the development of a fully three-dimensional (3D) free-surface flow model.However, the required computational fluid dynamics (CFD) simulation is time consuming and its validation requires experimental parameters that are difficult to determine.Since the goal was to obtain insights into the interaction of the droplet and the metallic target even after the onset of erosion caused by a relatively large number of impacts, other numerical methods that are faster and easier to implement were used in this work.Namely, the impact of individual droplets was studied with a combination of Smoothed-Particle Hydrodynamics (SPH) and Finite Element Method (FEM) modeling.The dynamics of the impacting water droplet were simulated within the SPH method, whereas the mechanical behavior of the metallic target was represented within the FEM framework.
Verma et al. recently applied this type of hybrid approach in a comprehensive simulation of the effect of rain droplets impacting the leading edge of wind turbine edges [28].The blades considered are made from a fiber composite protected with a coating of an epoxybased thermosetting polymeric resin.Diverse parametric studies, varying droplet sizes, impact angles, and impact velocities were considered.Peak impact forces and maximum effective plastic strain were calculated.However, due to the small finite element size and high number of SPH particles representing the rain droplets, among other complexities of the model, the simulations were limited to only ten impacts on the same location.Therefore, the onset of erosion damage on the coating was predicted analytically by performing an extrapolation based on known properties of the specific material [28].The approach to this work is different: Simple material/damage models together with relatively coarse and efficient finite element meshes and smoothed-particle method configurations were chosen to be able to simulate up to 1000 impacts on random positions on a representative area of the metallic surface.The novelty of the model lies in the capability of simulating, using reasonable computing resources, the complete incubation period and the onset of the water droplet erosion regime for certain ranges of impact angles and velocities.

Water Droplet Erosion Mechanisms
During impact, the droplet initially exerts normal pressure on the surface, which creates stress waves in a bulk material.S. S. Cook introduced the concept of "water hammer pressure" to describe the continuous pressure generated by droplet impacts on a solid surface [29].The water hammer pressure can be expressed as: where ρ 0 is the liquid density, C 0 is the speed of sound in the liquid, and V is the impact speed of the droplet.Field et al. modified this equation to take into account the shock wave properties of the solid surface [30], as the following equation shows: where ρ s and C s represent the density and speed of sound of the solid (target surface).
Heymann postulated the maximum peak value of the impact pressure [1]: which can be used as an estimation of the loading condition of each impact.If the generated stress wave inside the material has an amplitude high enough to overcome the dynamics fracture strength at certain locations, damage and fracture eventually can take place [11].
After releasing the impact pressure on the surface, the liquid droplet is expelled in all outward directions tangent to the surface, which is known as "water jetting", as shown in Figure 1a.On a perfectly flat surface, the water jetting effect initially does not cause erosion.However, above a certain impact velocity threshold, repeated impacts cause plastic deformations in the form of depressions (dimples) and the corresponding protruding edges on the surface (Figure 1b).As a result, water jetting interacts with these asperities, plastically shearing them and initiating micro-cracks on the surface [31].
plexities of the model, the simulations were limited to only ten impacts on the same location.Therefore, the onset of erosion damage on the coating was predicted analytically by performing an extrapolation based on known properties of the specific material [28].
The approach to this work is different: Simple material/damage models together with relatively coarse and efficient finite element meshes and smoothed-particle method configurations were chosen to be able to simulate up to 1000 impacts on random positions on a representative area of the metallic surface.The novelty of the model lies in the capability of simulating, using reasonable computing resources, the complete incubation period and the onset of the water droplet erosion regime for certain ranges of impact angles and velocities.

Water Droplet Erosion Mechanisms
During impact, the droplet initially exerts normal pressure on the surface, which creates stress waves in a bulk material.S. S. Cook introduced the concept of "water hammer pressure" to describe the continuous pressure generated by droplet impacts on a solid surface [29].The water hammer pressure can be expressed as: where  is the liquid density,  is the speed of sound in the liquid, and  is the impact speed of the droplet.Field et al. modified this equation to take into account the shock wave properties of the solid surface [30], as the following equation shows: where  and  represent the density and speed of sound of the solid (target surface).
Heymann postulated the maximum peak value of the impact pressure [1]: which can be used as an estimation of the loading condition of each impact.If the generated stress wave inside the material has an amplitude high enough to overcome the dynamics fracture strength at certain locations, damage and fracture eventually can take place [11].
After releasing the impact pressure on the surface, the liquid droplet is expelled in all outward directions tangent to the surface, which is known as "water je ing", as shown in Figure 1a.On a perfectly flat surface, the water je ing effect initially does not cause erosion.However, above a certain impact velocity threshold, repeated impacts cause plastic deformations in the form of depressions (dimples) and the corresponding protruding edges on the surface (Figure 1b).As a result, water je ing interacts with these asperities, plastically shearing them and initiating micro-cracks on the surface [31].The water hammer pressure component of subsequent droplet impacts produces more deformations and cracks on the already roughened surface and then damage (in the form of cracks, hardening, or/and fatigue) accumulates progressively.This roughened surface in turn alters the conditions for the water jetting and instead of a uniform spreading splashing of the droplet takes place (Figure 1c) [32].The splashing can be described as a crown formed by perturbations of the water film rim spreading on diverse angles reducing the amount and intensity of the water jetting effect.Therefore, there is a complex interplay among the water hammering, water jetting, and surface roughness effects that is necessary to be taken into account (Figure 1d) [33][34][35].This process of damage advancement in turn leads to crack growth, crack coalescence, fragmentation, pit formation, and measurable target mass removal (the end of the incubation period).

Modeling Methods
The interaction of the water droplets with the eroding surface has been simulated in this work using a hybrid Smoothed-Particle Hydrodynamics (SPH)/Finite Element Method (FEM) modeling scheme for multiple damage accumulation scenarios in detail.Namely, every individual droplet is represented by a number of virtual SPH particles impacting the target surface, resulting in the water hammer pressure during impact as well as the water jetting and its effects on the surrounding area [36].On the other hand, a target surface of Ti-6Al-4V titanium alloy was represented by a finite element mesh using the basic Johnson-Cook material constitutive model for both strength and failure [37][38][39].The Johnson-Cook model is broadly used due to its simplicity, relative numerical robustness, and availability of experimental values for empirical constants for most used metallic alloys.
The solid target material selected for the simulation was a titanium alloy known as Ti-6Al-4V.This annealed alloy has applications that range from wrought, cast, or additivemanufactured implants and prostheses, to parts and prototypes for the racing and aerospace sector, marine applications, the chemical industry, and gas turbines.The composition specifications and properties of Ti6Al4V are outlined in Tables 1 and 2: Table 1.The Ti6Al4V composition according to ASTM standards [40].

Johnson-Cook Material Constitutive Model for the Target
The Johnson-Cook (J-C) constitutive material model is usually employed in high dynamic process simulation and considers three effect factors: strain hardening, strain rate (viscosity), and thermal softening.In this model, the Mises flow stress, σ M , is calculated as shown in Equation ( 4).The first term in brackets represents the strain hardening factor.The expressions in the second and third terms in brackets represent the dynamic effects, which are functions of the plastic strain rate and the temperature, respectively, as: ε 0 are the equivalent plastic strain, the strain rate, and the reference strain rate, respectively.T, T 0 , and T f represent the temperature ( • C), the room temperature, and the melting temperature of the material.A, B, C, m, and n are constants, respectively, representing the static yield strength, the strain hardening modulus, the coefficient of the strain rate, the thermal-softening exponent, and the strain-hardening exponent.These parameters are determined based on the flow stress data obtained from mechanical tests for the specific material and temperature ranges [41], as shown in Table 3.

Johnson-Cook Failure Model for the Target
The Johnson-Cook damage criterion is based on the magnitude of the equivalent plastic strain and has been used to simulate material removal during the erosion process in this work.When the equivalent plastic strain reaches a predefined failure strain value, damage occurs and elements are removed.Equation ( 5) represents the equivalent plastic strain ε pl, f at material failure: where σ p is the pressure stress, σ M is the von Mises stress, and .
ε pl is the equivalent plastic strain rate.The constants d 1 , d 2 , d 3 , d 4 , and d 5 are material constants obtained from ElTobgy et al. [41] and listed in Table 4.According to [41].
The failure criterion is based on the equivalent plastic strain value, ε pl .Failure and removal on a specific finite element is assumed to occur when the damage variable D reaches or exceeds a value of 1, and is defined by the following expression: where ∆ε pl is an increment of the equivalent plastic strain and ε pl, f is the equivalent plastic strain at failure (Equation ( 5)).The summation process involves computing the total plastic strain accumulated over all increments during the analysis.

Smoothed-Particle Hydrodynamics (SPH) Method for Water Droplets
Smoothed-Particle Hydrodynamics (SPH) is a mesh-free Lagrangian method in which virtual particles without fixed connectivity carry material properties, which is used to simulate extreme deformation problems [42].The SPH method can be used in a wide range of applications such as quasi-incompressible fluid flow, viscous fluid flow, the high-velocity impact of solids, geophysical flows, fluid sloshing, wave engineering, ballistics, spraying, gas flow, and obliteration and fragmentation followed by secondary impacts [43,44].
The SPH method involves representing field variables such as density, velocity, and energy through partial differential equations.There are two key steps in its formulation: The kernel approximation and the particle approximation.In the kernel approximation, a function f (x) and its derivatives are expressed through an integral representation using a smoothing kernel function: where < f (x) > is the kernel approximation of the function f (x), x is the current position vector, x are the position vectors in the integral, Ω is the support domain, h is the smoothing length, and W is the smoothing function.The particle approximation involves discretizing the computational domain with a set of virtual particles that represent the fluid and its dynamical behavior.The continuous integral expression of Equation ( 7) is converted to a discretized form of all particles in the domain.Field variables are then approximated by summing the values of nearest neighbor particles.The particle approximation of the function at particle j is defined as: where N is the total number of particles in the region of influence of the current particle, m j is the mass of neighbor particle j, and ρ j is the density.
In order to apply the SPH to represent the dynamics of the water droplets impacting the metallic target, an "equation of state" for water is introduced.Such an equation of state is a semi-empirical relationship between pressure, volume, and temperature in the fluid and is approximately compatible with the conservation laws of hydrodynamics.In this work, the Mie-Grüneisen equation of state in its linear U s − U p Hugoniot form for the pressure p is obtained in terms of the internal energy per unit mass, where ρ 0 and c 0 are the density and the speed of sound in water, respectively, η is the nominal volumetric compressive strain, Γ 0 is a material constant, and s is the Hugoniot slope coefficient.For this work, a basic parameter set was chosen (see Table 5).The target material FEM model was a 3-dimensional solid deformable part with dimensions of 28 × 28 × 10 mm.The finite element type applied is the C3D8R brick element, which means a continuum (C), three-dimensional (3D), and eight-node reduced (8R) integration element.In order to obtain a realistic simulation of the erosion process on a representative surface, the target plate was divided into two parts: First, the "inner" part (dimensions 18 × 18 × 5 mm) on which the water droplet impacts effectively take place.Second, the surrounding "outer" part was used for reproducing realistic boundary conditions for the "inner" part (See Figure 2a).Namely, the lateral and bottom sides of the external outer part are maintained fixed ("encastre" boundary condition) to represent the fact that the simulated portion is part of a whole solid.The "outer" part is also big enough to absorb shock waves and deformations traveling from the "inner" box.It is important to reiterate that the outer part does not receive impacts and hence no element removal/erosion occurs there.
Metals 2023, 13, x FOR PEER REVIEW 7 of 13 gration element.In order to obtain a realistic simulation of the erosion process on a representative surface, the target plate was divided into two parts: First, the "inner" part (dimensions 18 × 18 × 5 mm) on which the water droplet impacts effectively take place.Second, the surrounding "outer" part was used for reproducing realistic boundary conditions for the "inner" part (See Figure 2a).Namely, the lateral and bo om sides of the external outer part are maintained fixed ("encastre" boundary condition) to represent the fact that the simulated portion is part of a whole solid.The "outer" part is also big enough to absorb shock waves and deformations traveling from the "inner" box.It is important to reiterate that the outer part does not receive impacts and hence no element removal/erosion occurs there.The exterior walls of the inner part and the interior walls of the "outer part" are tied together using a "tie constraint", which makes the translational and rotational motion equal for the two surfaces.Moreover, the "tie constraint" allows for dissimilar meshes for both parts as shown in Figure 2b.Therefore, the size of the finite elements on the "outer" part can be conveniently selected to be much larger (approximately 60 times on average) than the elements of the inner part.In addition, the mesh size distribution can be biased towards the vicinity of the boundary between both parts, saving many computational resources without compromising the accuracy of the calculation.
The size of the finite elements of the "inner part" is maintained small (edge size of the cubic finite element = 0.25 mm) in order to obtain convergence of the numerical computations.The inner part is composed of 103,600 elements whereas the outer part has only 6660 elements.Test simulations with even smaller finite elements were performed (edge size = 0.0125 mm) without noticeable qualitative differences in the results, indicating that convergence was granted.However, smaller finite elements always create problems with the available computational resources (computational time, storage, and analysis time) considering that mesh refinement should also be performed together with an increase in the number of SPH particles representing the water droplets.

Smoothed-Particle Hydrodynamics Simulation of the Water Droplets
To represent each individual water droplet within Abaqus FEA, a spherical geometry with a diameter of 4 mm was created composed initially from a special type of finite element named "continuum particle elements" PC3D (Figure A1a of Appendix A depicts 5 of those modeled water droplets).The relatively large diameter of the water droplet (4 mm) was chosen in order to obtain an incubation period within the computational resources available.However, for the cases of erosion of fan blades of aero-engines, the droplets have diameters within the 1-5 mm range [7].Each node of each PC3D element acts as an SPH virtual particle that interacts with all other SPH particles within a dominion The exterior walls of the inner part and the interior walls of the "outer part" are tied together using a "tie constraint", which makes the translational and rotational motion equal for the two surfaces.Moreover, the "tie constraint" allows for dissimilar meshes for both parts as shown in Figure 2b.Therefore, the size of the finite elements on the "outer" part can be conveniently selected to be much larger (approximately 60 times on average) than the elements of the inner part.In addition, the mesh size distribution can be biased towards the vicinity of the boundary between both parts, saving many computational resources without compromising the accuracy of the calculation.
The size of the finite elements of the "inner part" is maintained small (edge size of the cubic finite element = 0.25 mm) in order to obtain convergence of the numerical computations.The inner part is composed of 103,600 elements whereas the outer part has only 6660 elements.Test simulations with even smaller finite elements were performed (edge size = 0.0125 mm) without noticeable qualitative differences in the results, indicating that convergence was granted.However, smaller finite elements always create problems with the available computational resources (computational time, storage, and analysis time) considering that mesh refinement should also be performed together with an increase in the number of SPH particles representing the water droplets.

Smoothed-Particle Hydrodynamics Simulation of the Water Droplets
To represent each individual water droplet within Abaqus FEA, a spherical geometry with a diameter of 4 mm was created composed initially from a special type of finite element named "continuum particle elements" PC3D (Figure A1a of Appendix A depicts 5 of those modeled water droplets).The relatively large diameter of the water droplet (4 mm) was chosen in order to obtain an incubation period within the computational resources available.However, for the cases of erosion of fan blades of aero-engines, the droplets have diameters within the 1-5 mm range [7].Each node of each PC3D element acts as an SPH virtual particle that interacts with all other SPH particles within a dominion of influence whose radius is referred to as the smoothing length h, as defined in Section 3.3, Equations ( 7) and (8).During the numerical calculation, the SPH particles interact among themselves, mimicking a hydrodynamic behavior.In addition, the SPH particles also have a mechanical interaction with the finite elements representing the target surface by means of a "general contact" condition.This allows for the calculation of the interaction among the nodes of the C3D8R elements of the target and those from the SPH particles of the water droplet.Eventually, after a number of impacts, some finite elements of the target accumulate enough damage to reach the failure criterion (D = 1, see Equation ( 6)) and are removed from the surface exposing previously underlying elements to subsequent water droplet impacts, damage, and eventual removal.
The simulation of multiple impacts on a representative surface requires an algorithm that generates random impacts for a range of diverse impact angles and velocities (angles 10 • -90 • ).Appendix A presents a detailed discussion of the practical considerations for performing simulations of 1000 random impacts on a sequence.

Results and Discussion
Figure 3 summarizes the results of the simulation of impacting 1000 water droplets at angles ranging from 10 • to 90 • degrees (90 • being on the normal direction of the initial flat surface) at a velocity of 600 m/s on the area of the target defined in Appendix A (Figure A1b). Figure 3a shows that element removal, measured as volume loss in mm 3 , occurs only for angles above 60 degrees.At that specific angle, element removal occurs only after roughly 650 impacts, which could be considered the end of the incubation regime (designated as the "eir" variable for each angle in the figure).The erosion is much more effective for an impacting angle of 70 • .At angles of 80 • and 90 • , the onset and evolution of the erosion regime show similar features.Figure 3b shows the total volume loss after 1000 impacts at a velocity of 600 m/s for the whole range of impact angles.
of influence whose radius is referred to as the smoothing length ℎ, as defined in Section 3.3, Equations ( 7) and (8).During the numerical calculation, the SPH particles interact among themselves, mimicking a hydrodynamic behavior.In addition, the SPH particles also have a mechanical interaction with the finite elements representing the target surface by means of a "general contact" condition.This allows for the calculation of the interaction among the nodes of the C3D8R elements of the target and those from the SPH particles of the water droplet.Eventually, after a number of impacts, some finite elements of the target accumulate enough damage to reach the failure criterion ( = 1, see Equation ( 6)) and are removed from the surface exposing previously underlying elements to subsequent water droplet impacts, damage, and eventual removal.
The simulation of multiple impacts on a representative surface requires an algorithm that generates random impacts for a range of diverse impact angles and velocities (angles 10°-90°).Appendix A presents a detailed discussion of the practical considerations for performing simulations of 1000 random impacts on a sequence.

Results and Discussion
Figure 3 summarizes the results of the simulation of impacting 1000 water droplets at angles ranging from 10° to 90° degrees (90° being on the normal direction of the initial flat surface) at a velocity of 600 m/s on the area of the target defined in Appendix A (Figure A1b). Figure 3a shows that element removal, measured as volume loss in mm 3 , occurs only for angles above 60 degrees.At that specific angle, element removal occurs only after roughly 650 impacts, which could be considered the end of the incubation regime (designated as the "eir" variable for each angle in the figure).The erosion is much more effective for an impacting angle of 70°.At angles of 80° and 90°, the onset and evolution of the erosion regime show similar features.Figure 3b shows the total volume loss after 1000 impacts at a velocity of 600 m/s for the whole range of impact angles.Figure 4a shows how at an impact velocity of 800 m/s the onset of erosion occurs at impact angles larger than 30° and the number of impacts needed for the onset of erosion are lower than in the previous case.Figure 4b shows a maximum of erosion at 80°. Figure 4a shows how at an impact velocity of 800 m/s the onset of erosion occurs at impact angles larger than 30 • and the number of impacts needed for the onset of erosion are lower than in the previous case.Figure 4b shows a maximum of erosion at 80 • .
As expected, at an impact velocity of 1000 m/s, the incubation regime is much shorter and starts at a low impact angle of 30 • (Figure 5a).The maximum volume loss occurs at an impact angle of 70 • degrees (Figure 5b).The impact velocity v =1000 m/s implies an impact kinetic energy 2.777 times larger than the corresponding energy at 600 m/s (due to the kinetic energy definition: K = mv 2 /2).However, the volume loss does not grow linearly with the velocity nor the kinetic energy for the selected range of velocities.Namely, for v = 600 m/s, the volume loss at 80 • after 1000 impacts is 21.150 mm 3 , whereas the corresponding volume loss for v = 1000 m/s is 211.172mm 3 (approximately 10 times larger).This is due to the strong non-linear dynamics of the erosion process and the complex mechanisms that emerge when impact velocity and angle are varied.As expected, at an impact velocity of 1000 m/s, the incubation regime is much shorter and starts at a low impact angle of 30° (Figure 5a).The maximum volume loss occurs at an impact angle of 70° degrees (Figure 5b).The impact velocity  =1000 m/s implies an impact kinetic energy 2.777 times larger than the corresponding energy at 600 m/s (due to the kinetic energy definition:  =  /2).However, the volume loss does not grow linearly with the velocity nor the kinetic energy for the selected range of velocities.Namely, for  = 600 m/s, the volume loss at 80° after 1000 impacts is 21.150 mm 3 , whereas the corresponding volume loss for v = 1000 m/s is 211.172mm 3 (approximately 10 times larger).This is due to the strong non-linear dynamics of the erosion process and the complex mechanisms that emerge when impact velocity and angle are varied.There is an important detail to be considered in Figure 5a,b: The simulation for an impact angle of 90° stopped shortly after 900 impacts due to numerical instability.There-   3 versus the accumulated number of water droplet impacts for impact velocity  = 800 m/s at different impact angles.The variable "eir" (from "end of incubation regime") is approximately the number of impacts required for erosion for each impact angle.For angles lower than 40°, there is no erosion.(b) The total volume loss after 1000 impacts for different impact angles.See a detailed description in the text.
As expected, at an impact velocity of 1000 m/s, the incubation regime is much shorter and starts at a low impact angle of 30° (Figure 5a).The maximum volume loss occurs at an impact angle of 70° degrees (Figure 5b).The impact velocity  =1000 m/s implies an impact kinetic energy 2.777 times larger than the corresponding energy at 600 m/s (due to the kinetic energy definition:  =  /2).However, the volume loss does not grow linearly with the velocity nor the kinetic energy for the selected range of velocities.Namely, for  = 600 m/s, the volume loss at 80° after 1000 impacts is 21.150 mm 3 , whereas the corresponding volume loss for v = 1000 m/s is 211.172mm 3 (approximately 10 times larger).This is due to the strong non-linear dynamics of the erosion process and the complex mechanisms that emerge when impact velocity and angle are varied.There is an important detail to be considered in Figure 5a,b: The simulation for an impact angle of 90° stopped shortly after 900 impacts due to numerical instability.Therefore, the red curve with circle symbols in Figure 5a does not reach the "number of impacts" limit of 1000.Correspondingly, the bar at 90° in Figure 5b is missing.This numerical instability is due to the high velocity of 1000 m/s in the vertical direction and the re- There is an important detail to be considered in Figure 5a,b: The simulation for an impact angle of 90 • stopped shortly after 900 impacts due to numerical instability.Therefore, the red curve with circle symbols in Figure 5a does not reach the "number of impacts" limit of 1000.Correspondingly, the bar at 90 • in Figure 5b is missing.This numerical instability is due to the high velocity of 1000 m/s in the vertical direction and the resulting extreme deformations, which cannot be solved numerically using the explicit algorithm.For future work, the most obvious solution would be to refine the target mesh and increase the number of SPH particles representing the water droplets.As mentioned before, this would require much higher computational resources.

Conclusions
In this study, we used a well-established software for finite element analysis (Abaqus FEA) to simulate water droplet erosion on a Ti6Al4V target by coupling the Smoothed-Particle Hydrodynamics (SPH) and the Finite Element Method (FEM) approaches.The main achievement is to have combined simple material and damage models to investigate the target surface and particle-based dynamical behavior of the water droplets during impacts, which enabled the simulation of the incubation regime.Moreover, the designed simulation box (with "inner" and "outer" parts) and the optimized mesh geometry allowed us to produce a physically sound representation of the boundary conditions required for exploring the effect of a wide range of impact velocities and angles.The mesh size significantly influenced the accuracy of the simulation results, as usual in finite element analyses.As explained in the previous section, a finite element size of 0.25 mm was selected ("inner" part) as an optimal compromise between simulation accuracy and the availability of computational resources.
The simulation results shown in Figures 3-5 reproduce the main qualitative aspects of water droplet erosion of metallic surfaces reported in the literature [3,11,13]:

•
For impact velocities below certain threshold, erosion eventually ceases to occur (within a reasonable lifetime of the metallic targets and independent of any impact angle).

•
Above the threshold velocity mentioned before, an incubation regime exists where neither mass removal nor surface roughening occurs.The duration of the incubation regime and the onset of erosion (in terms of the number of accumulated impacts) decreases with increasing impact velocity/kinetic energy.

•
Depending on the impact velocity, there is an impact angle lower than 90 • where the volume loss is maximum.Namely, in our simulations, the maximum is 80 • for velocities v = 600 m/s and 800 m/s and 70 • for 1000 m/s • In spite of these achievements, it is important to emphasize that the model used in this work is a simplification of the physical mechanisms taking place during real water droplet erosion.Moreover, the specific ranges of simulation parameters selected in this work, namely large water droplet sizes (4 mm diameter) and high velocities (600-1000 m/s), allowed us to obtain erosion using reasonable computational resources.However, for this very reason, it is difficult to conduct a direct quantitative comparison with experiments because there are currently no data available on the duration of the incubation period for the combination of ranges of droplet sizes and impact velocities used in this study (to the best of our knowledge).

•
Further work should focus on applying more sophisticated material and damage models, taking into account the influence of temperature and corrosion among other variables.Moreover, more realistic combinations of droplet size and impact velocity ranges should be used in order to compare with the experimental results of the literature.A second goal would be to decrease the size of the finite elements and increase the number of SPH particles representing the water droplet.Finally, it is worth considering the application of an SPH model (with the relevant material and damage model) to the impacted metallic surface itself in order to produce a more realistic evolution of the extreme deformations that occur on the real surfaces of material alloys.In order to simulate 1000 impacts, 200 simulations of the impact of five particles random positions were performed in sequence (see Figure A2).Each one of the 200 sim lations first loads the finite element mesh of the target computed on the previous simu tion (which includes accumulated plastic deformation, stress damage, and removed e ments) and then the impact of five new droplets on random positions (within the r square area described before) is performed.This allows for efficient use of the availa data storage and computer memory and eases tracking of the evolution of the dama and number of removed elements during the simulation of 1000 impacts.In order to simulate 1000 impacts, 200 simulations of the impact of five particles on random positions were performed in sequence (see Figure A2).Each one of the 200 simulations first loads the finite element mesh of the target computed on the previous simulation (which includes accumulated plastic deformation, stress damage, and removed elements) and then the impact of five new droplets on random positions (within the red square area described before) is performed.This allows for efficient use of the available data storage and computer memory and eases tracking of the evolution of the damage and number of removed elements during the simulation of 1000 impacts.

Figure 1 .
Figure 1.A schematic representation of the physical processes modeled.(a) The droplet impact on an initially flat and intact surface produces "water hammer" pressure and water je ing (in perpendicular direction as indicated by the vertical blue arrow).(b) After a number of impacts, internal stresses, damage, and plastic deformation lead to dimple formation.(c) Repeated water je ing shears the rims of the dimples and their related asperities, which adds to the accumulated damage

Figure 1 .
Figure 1.A schematic representation of the physical processes modeled.(a) The droplet impact on an initially flat and intact surface produces "water hammer" pressure and water jetting (in perpendicular direction as indicated by the vertical blue arrow).(b) After a number of impacts, internal stresses, damage, and plastic deformation lead to dimple formation.(c) Repeated water jetting shears the rims of the dimples and their related asperities, which adds to the accumulated damage and fracture.(d) At the advanced stage of the incubation period, the surface is roughened, thus altering the effectivity of the droplet impacts.See a detailed description in the text.

Figure 2 .
Figure 2. (a) The finite element simulation box geometry of the metallic target.(b) The target is composed of the "inner" part (green box) and "outer" part (gray volume) with the shown dimensions.The inner and outer parts are tied together by means of a "tie" constraint.The simulated droplets only impact the "inner" part.See a detailed description in the text.

Figure 2 .
Figure 2. (a) The finite element simulation box geometry of the metallic target.(b) The target is composed of the "inner" part (green box) and "outer" part (gray volume) with the shown dimensions.The inner and outer parts are tied together by means of a "tie" constraint.The simulated droplets only impact the "inner" part.See a detailed description in the text.

Figure 3 .
Figure 3. (a) The volume loss measured in mm 3 versus the accumulated number of 4 mm diameter water droplet impacts for impact velocity  = 600 m/s at different impact angles (for angles lower than 60° there is no erosion.The variable "eir" (from "end of incubation regime") is approximately the number of impacts required for the erosion for each impact angle.(b) The total volume loss after 1000 impacts for different impact angles.See a detailed description in the text.

Figure 3 .
Figure 3. (a) The volume loss measured in mm 3 versus the accumulated number of 4 mm diameter water droplet impacts for impact velocity v = 600 m/s at different impact angles (for angles lower than 60 • there is no erosion.The variable "eir" (from "end of incubation regime") is approximately the number of impacts required for the erosion for each impact angle.(b) The total volume loss after 1000 impacts for different impact angles.See a detailed description in the text.

Figure 4 .
Figure 4. (a) The volume loss measured in mm3 versus the accumulated number of water droplet impacts for impact velocity  = 800 m/s at different impact angles.The variable "eir" (from "end of incubation regime") is approximately the number of impacts required for erosion for each impact angle.For angles lower than 40°, there is no erosion.(b) The total volume loss after 1000 impacts for different impact angles.See a detailed description in the text.

Figure 5 .
Figure 5. (a) The volume loss measured in mm 3 versus the accumulated number of water droplet impacts for impact velocity  = 1000 m/s at different impact angles.The variable "eir" (from "end of incubation regime") is approximately the number of impacts required for erosion for each impact angle.For angles lower than 30°, there is no erosion.(b) The total volume loss after 1000 impacts for different impact angles.See a detailed description in the text.

Figure 4 .
Figure 4. (a) The volume loss measured in mm 3 versus the accumulated number of water droplet impacts for impact velocity v = 800 m/s at different impact angles.The variable "eir" (from "end of incubation regime") is approximately the number of impacts required for erosion for each impact angle.For angles lower than 40 • , there is no erosion.(b) The total volume loss after 1000 impacts for different impact angles.See a detailed description in the text.

Figure 4 .
Figure 4. (a) The volume loss measured in mm3 versus the accumulated number of water droplet impacts for impact velocity  = 800 m/s at different impact angles.The variable "eir" (from "end of incubation regime") is approximately the number of impacts required for erosion for each impact angle.For angles lower than 40°, there is no erosion.(b) The total volume loss after 1000 impacts for different impact angles.See a detailed description in the text.

Figure 5 .
Figure 5. (a) The volume loss measured in mm3 versus the accumulated number of water droplet impacts for impact velocity  = 1000 m/s at different impact angles.The variable "eir" (from "end of incubation regime") is approximately the number of impacts required for erosion for each impact angle.For angles lower than 30°, there is no erosion.(b) The total volume loss after 1000 impacts for different impact angles.See a detailed description in the text.

Figure 5 .
Figure 5. (a)The volume loss measured in mm 3 versus the accumulated number of water droplet impacts for impact velocity v = 1000 m/s at different impact angles.The variable "eir" (from "end of incubation regime") is approximately the number of impacts required for erosion for each impact angle.For angles lower than 30 • , there is no erosion.(b) The total volume loss after 1000 impacts for different impact angles.See a detailed description in the text.

Figure A1 .
Figure A1.(a) Five water droplets (represented by a spherical volume of SPH particles with a 4 m diameter) almost simultaneously impact the FEM simulation box representing the metallic tar surface.(b) The top view of the xy-plane.In order to avoid damage or removal of finite elements the "outer" part, random impacts occur only on a smaller area of the "inner" part shown by a square.Magenta horizontal lines show the eventual maximum range of eventual damage/elem removal.(c) A view of the xz-plane of the same five impacting particles.The dashed lines indic the direction of the impacts (impact angle of 30°).See a detailed description in the text.

Figure A1 .
Figure A1.(a) Five water droplets (represented by a spherical volume of SPH particles with a 4 mm diameter) almost simultaneously impact the FEM simulation box representing the metallic target surface.(b) The view of the xy-plane.In order to avoid damage or removal of finite elements of the "outer" part, random impacts occur only on a smaller area of the "inner" part shown by a red square.Magenta horizontal lines show the eventual maximum range of eventual damage/element removal.(c) A view of the xz-plane of the same five impacting particles.The dashed lines indicate the direction of the impacts (impact angle of 30 • ).See a detailed description in the text.

Figure A2 .
Figure A2.Three temporal frames of an individual SPH/FEM simulation of the impact of five water droplets at the same velocity of 1000 m/s and impact angle of 70 • on random positions as described in FigureA1.Please note that only the "inner part" of the target is being depicted here (dimensions 18 × 18 × 5 mm, see Figure2a).(a) At t = 0 s, the simulation starts loading the target with the deformations, stress, damage, and removed elements of previous impacts and a new set of five droplets moving toward the surface.(b) At t = 2.9 ms, the five droplets are interacting with the eroded target.(c) Finally, at t = 8.7 ms, the SPH particles representing the droplet rebound off the surface.The simulation of the next five droplets starts with this final state target FEM mesh.

Table 3 .
The Johnson-Cook constitutive material model parameters of the Ti6A4V alloy considered in this work.

Table 4 .
The Johnson-Cook failure model parameters of Equation (5) used for the simulations.

Table 5 .
The Mie-Grüneisen equation of state parameters of Equation (9) used in this work for the SPH modeling of impacting droplets.
4. Numerical Simulations of Multiple Droplet Impacts4.1.Finite Element Numerical Simulation of the Target Surface Metals 2023, 13, x FOR PEER REVIEW 9 of 13