Shrey Joshi
email: shrey.joshi@simap.grenoble-inp.fr

Jean Pierre Franc

Giovanni Ghigliotti

Marc Fivel

An axisymmetric Solid SPH solver with consistent treatment of particles close to the symmetry axis

Keywords: Smoothed Particle Hydrodynamics, plasticity, cavitation, material damage, axisymmetric

A meshless Smoothed Particle Hydrodynamics solid solver is developed in order to study fluidstructure interactions and to predict cavitation erosion. The solid solver is developed in-house in an axisymmetric configuration. The existing SPH methods dedicated to solid materials do not allow a consistent treatment of particles close to the symmetry axis, so a density correction scheme is proposed here to derive new density and momentum equations for solid mechanics in axisymmetric SPH formulation. The new SPH equations are coded in the solver and the SPH solid solver is then validated against FEM results which shows excellent agreement.

Introduction

Cavitation erosion is a major issue amongst a wide range of equipment's such as hydraulic devices, diesel injectors, artificial heart valves, etc.. Cavitation is defined as the appearance of vapor cavities inside a continuous and homogeneous liquid medium. The generation of vapor cavities could happen due to various reasons, but the cause is mostly associated to the drop of local pressure below the vapor pressure. Figure 1 shows a pressure-temperature phase diagram of water, showing that a drop of pressure below the vapor pressure leads to vaporization of liquid. When the vapor cavities collapse is located near the surface, it can lead to the formation of a high velocity micro jet and subsequent shock wave. Such an asymmetry in the collapse dynamics is due to a weaker fluid flow on the wall side (lower flow velocity due to the presence of the wall). Figure 2 shows a schematic description of a collapsing bubble and the resulting micro-jet. During the bubble collapse, the solid surface experiences a high pressure due to the shock wave and the high velocity micro jet hitting the surface. The magnitude of the pressure acting on the surface could be high enough to cause plastic damage to the material. Although the experimental measurement of these pressures acting on the surface still remains a challenge, various studies have estimated it to be around a few GPa [START_REF] Philipp | Cavitation erosion by single laser-produced bubbles[END_REF][START_REF] Vogel | Optical and acoustic investigations of the dynamics of laser-produced cavitation bubbles near a solid boundary[END_REF][START_REF] Ahmed | Fatigue failure of SUS 304 caused by vibratory cavitation erosion[END_REF][START_REF] Jones | An experimental study of the forces generated by the collapse of transient cavities in water[END_REF][START_REF] Momma | A study of pressures and erosion produced by collapsing cavitation[END_REF]. Moreover, experimental investigations offer only limited information about bubble collapses and the subsequent formation of micro jets and pressure waves since any intrusive measurement can cause deviations. Such difficulties in experimental investigations have led to the use of numerical studies to understand and analyze cavitation.

Most numerical studies presented in the past few years have focused on solving the dynamics of the cavitation bubble collapse inside the fluid. However, more work need to be done in order to get a thorough understanding of material damage induced by collapsing bubbles. For instance, most numerical CFD studies estimate cavitation erosion-prone areas using either peak pressures or Cavitation Aggressiveness Index (CAI) [START_REF] Phoevos | A cavitation aggressiveness index within the Reynolds averaged Navier Stokes methodology for cavitating flows[END_REF][START_REF] Li | Assessment of cavitation erosion with a multiphase Reynolds-Averaged Navier-Stokes method[END_REF][START_REF] Wang | Numerical Computation of Shock Waves in a Spherical Cloud of Cavitation Bubbles[END_REF], however it is still debatable as to whether peak pressure or CAI can accurately characterize the cavitation erosion process taking place in the material. Hence to address the above issues better, recently there has been a thrust towards modelling cavitation erosion using Fluid-Structure Interaction solvers to understand the phenomenon of cavitation erosion in entirety [START_REF] Hsiao | Modelling of material pitting from cavitation bubble collapse[END_REF][START_REF] Pöhl | Evaluation of cavitation-induced pressure loads applied to material surfaces by finite-element-assisted pit analysis and numerical investigation of the elasto-plastic deformation of metallic materials[END_REF][START_REF] Turangan | Numerical studies of cavitation erosion on an elasticplastic material caused by shock-induced bubble collapse[END_REF]. A thorough numerical study of cavitation that can provide a holistic understanding, ideally requires a two-way fluid-structure interaction coupling to get realistic results for cavitation erosion. A comprehensive cavitation numerical model should ideally include the following elements in the solver: a fluid solver, a solid solver and a fluid structure interaction scheme.

Accepted Version of : Joshi, S., Franc, J.P., Ghigliotti, G. et al. An axisymmetric solid SPH solver with consistent treatment of particles close to the symmetry axis. Comp. Part. Mech. (2020). https://doi.org/10.1007/s40571-019-00310-8 *Corresponding Author, Shrey Joshi: shrey.joshi@simap.grenoble-inp.fr The problem is quite complex since the fluid solver (typically a Eulerian Finite Volume Method (FVM) code) and the solid solver (typically a Lagrangian Finite Element Method (FEM) code) are generally two different codes which then need to be coupled in order to communicate data. Moreover, the elasto-plastic deformation of the solid medium due to cavitation loads implies the fluid mesh has to be re-constructed to match the solid mesh and maintain the continuity at the interface. Conventionally, ALE (Arbitrary Lagrangian Eulerian) methods have been used to simulate such a problem [START_REF] Paquette | Interaction Fluide-Structure en Érosion de Cavitation[END_REF]. To overcome the above complexity, in the present study a first attempt has been made to solve cavitation erosion using a meshless particle method, namely the Smoothed Particle Hydrodynamics (SPH), to model both the fluid and the solid behavior in a unique Lagrangian framework. The method offers the following advantages:

• Both fluid and solid response can be captured using the same solver and the same numerical method i.e. SPH. Both the solid and the fluid can be coded and solved within a single code making it much easier from the development perspective.

• Coupling a FVM solver to a FEM solver can be quite complicated and requires dedicated procedures for data transfer across the two codes. The problem is eliminated here since the solver uses SPH for both solid and fluid in a common code; no such coupling is required; the data are easily exchanged within the RAM memory. For each particle, one only needs to specify if it is a fluid or a solid particle.

• The tricky mesh reconstruction for the fluid domain is eliminated when using meshless methods like SPH where both the fluid and the solid are solved using a Lagrangian formulation and hence such deformation are already taken care of via particle movement. To develop the SPH FSI solver, the existing 2D open source fluid SPH code SPHYSICS is used as a basis to develop the solver further [START_REF] Gómez-Gesteira | User Guide for the SPHysics Code v2.0[END_REF]. However, a 2D fluid structure interaction (FSI) solver for a collapsing bubble over a solid medium does not represent a spherical cavity but rather an infinitely long cylindrical cavity which is a misrepresentation of the geometry. Also the response of a 2D solid medium is quantitatively different from that in the actual case. To make sure the cavitation bubbles and the solid medium are represented as in reality, the 2D solver is modified to 2D axisymmetric. However, axisymmetric SPH solver suffer from an inconsistent definition of density near the symmetry axis. This paper presents a new method to mathematically resolve this issue inspired by an approach used in fluid axisymmetric SPH simulation. New density and momentum equations are derived for solid axisymmetric SPH solver which can treat particles to obtain accurate density field close to the axis, the solver is finally validated against FEM results.

SPH method a. SPH fundamentals

Smoothed Particle Hydrodynamics (SPH) is a meshless numerical method to solve partial differential equations by discretizing the computational domain using set of particles. It was first introduced by Lucy [START_REF] Lucy | A numerical approach to the testing of the fission hypothesis[END_REF] in 1977 for astrophysics calculations and later by Monaghan [START_REF] Monaghan | Artificial viscosity for particle methods[END_REF] for fluid flows. The method is based around interpolations, where a quantity can be estimated at any point in space using its values at a set of neighboring particles by using a kernel function.

The integral interpolant of any function A at a position r is defined by the integration over the entire domain (Ω):

  ' ) , ' ( ) ' ( dr h r r W r A A     r ( 1
)
where W is the kernel function and h an associated length. The kernel function (W) is used for interpolation and must verify the following properties: firstly, the integration of the kernel function over the whole domain should lead to unity as given by equation 2 and secondly, as the limit of h (commonly known as smoothing length) tends to zero the kernel function becomes a Dirac-Delta function as given by equation 3.

1 = h)dr' , r' - W(r  (2) and ) r' - (r = h) , r' - W(r lim 0   h (3)
Several formulations have been proposed for the Kernel function in the past and have been used in SPHYSICS [START_REF] Gómez-Gesteira | User Guide for the SPHysics Code v2.0[END_REF]. In the present paper we chose to use a cubic spline function as detailed in appendix.

In order to use equation 1 for numerical purpose, the equation is transformed into a discretized form where it takes the form of a summation interpolant given by,

      j j N j j j m h W A A      1 , r r r r (4)
where the summation index j denotes a particle label, and the summation is performed over all the particles. Particle j has mass mj, position rj and density ρj. The value of any quantity A at rj is denoted by A(rj).

The key element is that we can write a differentiable interpolant of a function from its values at the particles by using a kernel which is differentiable. Derivatives of this kernel can be obtained by ordinary differentiation; there is no need to use finite differences or any other numerical method. For instance, if we want to compute the gradient ∇A, we can use

      j j N j j j m h , r r W r A r A        1     (5)

b. 2D SPH formulation for Solid mechanics

Accepted Version of : Joshi, S., Franc, J.P., Ghigliotti, G. et The continuity and the momentum equations are formulated in SPH form by transforming the equations into a summation over a set of discrete particles, interpolated using the kernel function. The density equation takes the form

  b ab b a W m  (6)
where a  is the density of particle a, b is the index for neighboring particles, b m is the mass and ab W is the value of kernel function centered at particle a and estimated at the neighboring particle, b. Within the framework of isotropic linear elasticity, the stress state at a given time step for a given particle is calculated from the constitutive equation,

kk ij ij ij µ       2 (7)
where σ is the stress tensor, ε is the strain tensor,  and  the Lame elastic constants.

Once density and stress are calculated from the above equations, the values are used in the SPH form of the momentum equation given below,

22 i ij ij a a b ab b ab j b a b a dv dW m dt dx           (8)
where Пab denotes artificial viscosity and

1 1 ( ) ( ) 2 n ij i i j j ab ab a a b b a b a ji b aa dW dW m v v v v dx dx           (9)
For modelling strain-rate dependent plasticity, we use the phenomenological Johnson-Cook model [START_REF] Johnson | Fracture characteristics of three metals subjected to various strains, strain rates, temperature and pressure[END_REF] In this paper, the temperature effects are neglected and equation [START_REF] Pöhl | Evaluation of cavitation-induced pressure loads applied to material surfaces by finite-element-assisted pit analysis and numerical investigation of the elasto-plastic deformation of metallic materials[END_REF]

reduces to, )] ln( 1 ][ ) ( [ ) , , ( * p n p p p y C B A T           (11)
In the Johnson-Cook model, A, B, C, n and m are material constants. Also 0 p  is the reference plastic strain-rate of the quasi- static test used to determine the yield and hardening parameters A, B and n. T0 is a reference temperature, and Tm is the reference melting temperature. For conditions where T * < 0, we assume that m = 1. Details regarding plasticity calculation and return mapping algorithm for non-linear isotropic hardening is provided in Appendix B.

c. Moving the particles

The particles are moved using XSPH variant [START_REF] Monaghan | On the problem of penetration in particle methods[END_REF] Accepted Version of : Joshi, S., Franc, J.P., Ghigliotti, G. et . The method moves the particle with an average velocity (weighted average based on parameter ε) of the particles in the neighborhood, this helps in avoiding inter-penetration of particles.

d. Time integration

Consider the momentum equation (equation 13) and position equation (equation 14) in the following form

a a F dt v d    ( 13
) a a V dt r d    ( 14
)
where Fa represents the force and Va represents the velocity contribution from particle a and from neighboring particles (XSPH correction).

The predictor step uses the time derivatives from the last time step to predict velocities and position at half-time step

n a n a n a F t v v    2 2 / 1     (15) n a n a n a V t r r    2 2 / 1     (16)
These values are then corrected using derivatives estimated at the half step

2 / 1 2 / 1 2      n a n a n a F t v v    ( 17
) 2 / 1 2 / 1 2      n a n a n a V t r r    (18)
Finally using the values of velocity and position from the predictor and corrector step one can obtain the values at the end of the time step by,

n a n a n a v v v        2 / 1 1 2 (19) n a n a n a r r r        2 / 1 1 2 (20)

e. Artificial viscosity

The artificial viscosity term to be used in equation 8 is given by, additional information regarding the equations 21 and 22 can be found in the work by Monaghan [START_REF] Monaghan | Artificial viscosity for particle methods[END_REF] Accepted Version of : Joshi, S., Franc, J.P., Ghigliotti, G. et al. An axisymmetric solid SPH solver with consistent treatment of particles close to the symmetry axis. Comp. Part. Mech. (2020). https://doi.org/10.1007/s40571-019-00310-8 *Corresponding Author, Shrey Joshi: shrey.joshi@simap.grenoble-inp.fr

          0 0 0 ab ab ab ab ab ab ab ab r v r v µ c       (21) 2 2    ab ab ab ab r r v h µ    (22)
Please note that the symbol µ in equation 7 is different from the symbol µab in equation 21-22.

Axisymmetric SPH for Solid Mechanics

a. Methodology

The original solver SPHYSICS_2D only solves for fluid in 2D (3D SPHYSICS solver is also available however this study used the 2D version of the code), a solid solver is added to the original fluid SPHYSICS in 2D plane strain. However, in that case, a 2D bubble represents a 3D infinitely long cylindrical bubble which is far from reality. Also a 2D plane-strain solid simulation is physically very different from the reality (in 3D) even for a similar cavitation load. To overcome this and not increase the simulation time significantly we chose to convert the 2D solver to an axisymmetric solver so that the spherical shape of the bubble is preserved. We present a detailed methodology which can be adopted to change 2D SPH solid solver to axisymmetric. Broadly the following changes are required to change the solver to axisymmetric:

1. Changing formulations and equations to cylindrical coordinate system 2. Adding a symmetry axis to the solver

Density corrections for particles close to the symmetry axis

An elegant approach to the axisymmetric Euler fluid equations was proposed by Brookshaw [START_REF] Brookshaw | Smooth particle hydrodynamics in cylindrical coordinates[END_REF] who derived the SPH form of these basic equations using the minimal action principle (see work by Monaghan [START_REF] Monaghan | Smoothed Particle Hydrodynamics[END_REF] and references thereafter for the history of variational SPH). A similar approach is used in the present work and the following changes are made to convert the 2D Solid SPH solver to an axisymmetric Solid SPH solver (this can also be used as a guide to develop an axisymmetric solver from any commonly available 2D open source SPH solvers): the density in axisymmetric formulation can be classified as the 3D density (which is the actual material density) and a 2D effective density which is used in the axisymmetric calculations. The 2D effective density is given by,

   N b ab a a W m 1  (23)
where the 2D density is defined as

a a a r    2 
which of course is different from the 3D density a  . To achieve this the mass of particle is allocated according to the following equation,

drdz r m a a a   2  (24)
The z-axis is defined as the symmetry axis. As shown in Figure 3b particles at a distance of 2h or less from the axis are mirror imaged across the z-axis to create ghost boundary particles. For any particle i with coordinate z) (r, , velocity ) z , r (   , mass mi, density  i and stress . Thus position, velocity as well as other quantities such as density and stress for ghost particles are updated at each step not using the SPH equations but from the evolution of real particles. Even though the use of ghost particles is not strictly necessary in axisymmetric geometry, it is possible to use a line of reflective particle at the symmetry axis without using ghost particle. However, it is valuable to use ghost particles to correctly represent the density and its derivatives near the singularity axis.

The equations of motion in cylindrical coordinate system are given by, Accepted Version of : Joshi, S., Franc, J.P., Ghigliotti, G. et

                         1 1 1 rz rr rr (25) z z b r r z dt dv                 rz rz zz 1     (26)
where b is the body force vector.

The SPH form of the above equations of motion in cylindrical form are given by the following [START_REF] Wang | Stable axisymmetric SPH formulation with no axis singularity[END_REF],

              2 2 2 2 2 2 2                        (27)
                      2 2 2 2 2 2           (28)
To summarize, we point out the fundamental differences between a 2D/3D and an axisymmetric SPH method. Unlike in 2D or 3D SPH, in axisymmetric SPH the mass of the particle depends on the initial radius of the particle from the symmetry axis as given by equation 24. This is done so that each particle represents the mass of the entire ring (as show in figure 3a). A dummy density called 2D density is defined in axisymmetric SPH which represents the numerical density in the simulation. Since the 2D density does not change significantly during a fluid flow in weak compressible SPH, this in general means that the particle adjust their volume while they move towards or away from the symmetry axis. For e.g. if a particle in a flow starts to move towards the symmetry axis, the particles will automatically adjust the inter-particle distance such that the higher mass of the particle does not numerically induce a higher density values. This would imply that the particle volume (in 2D) in the simulation will increase while it moves towards the axis and vice versa. This however presents a challenge in defining the kernel smoothing length while the particles can adjust their interpartcile distances. To overcome that, the smoothing length is adjusted every few time steps for each particle such that the number of particles inside the kernel remains approximately the same as before. However, numerical issues near the symmetry axis are quite common since the standard kernel interpolation provides erroneous density values close to the symmetry axis. This can be solved using correction strategies near the symmetry axis as discussed in section 3b.

b. Correction close to symmetry axis

A major issue with axisymmetric SPH is the treatment of particles that are close to the symmetric axis. Indeed, large errors in density and consequently other quantities can be observed with particles at distance less than 2h from the axis. Such errors in density and other parameters such as velocity, stress, strain etc. close to the axis can lead inaccurate or even unstable simulations. A mathematically consistent solution to tackle this issue is especially necessary for a problem like cavitation where the dynamics of the imploding bubble, near the symmetry axis, is very rapid and therefore the region close to the symmetry axis is of utmost importance. To understand the issue better, let's consider a particle in the vicinity of the symmetry axis moving towards the axis. Equation 24 suggests that the 2D density of a particle is directly and linearly proportional to the radial distance from the axis. As the particle approaches the symmetry axis, the 2D density of the particle should tend to zero. However, because of the presence of the ghost particles, the profile of η = 2π|r|ρ, is no longer linear when the particle radius is lower than 2h. The limited capability of standard kernels to interpolate accurately non-linear functions introduces error in the density calculation. In other words, the density of a particle on the symmetry axis should be zero according to equation 24; however, the interpolation using particles on the right of the axis (actual particles) and the left of the axis (ghost particle) using the standard kernels would lead to a non-zero 2D density which is clearly wrong. Usually, the errors are small and the interpolation is precise to second order in h far from the axis. Unfortunately, close to the symmetry axis errors are much larger so that density and other physical quantities are not well reproduced.

One mathematically consistent way to overcome the issue is to use modified interpolation kernel according to the particular geometry of the system (spherical or cylindrical) [START_REF] Omang | SPH in spherical and cylindrical coordinates[END_REF]. The resulting scheme does not have inaccurate density when the Accepted Version of : Joshi, S., Franc, J.P., Ghigliotti, G. et al. An axisymmetric solid SPH solver with consistent treatment of particles close to the symmetry axis. Comp. Part. Mech. (2020). https://doi.org/10.1007/s40571-019-00310-8 *Corresponding Author, Shrey Joshi: shrey.joshi@simap.grenoble-inp.fr particles approach the axis. However, the modified kernel does not have an analytical expression and numerical integration is required to calculate the value of the kernel and its derivatives at every time step which require significantly higher computational effort. Later, a fitting formulae for the modified kernels was proposed [START_REF] Omang | Shock collisions in 3D using an axisymmetric regularized smoothed particle hydrodynamics code[END_REF] but still involving a large number of operations which slow down the calculation. An alternative approach was proposed where without modifying the basic SPH scheme, correction terms are introduced to equations 6, 25 and 26, which become significant only close to the symmetry axis [START_REF] García-Senz | Axisymmetric smoothed particle hydrodynamics with selfgravity[END_REF]. A detailed derivation has been provided by García-Senz et al. [START_REF] García-Senz | Axisymmetric smoothed particle hydrodynamics with selfgravity[END_REF] where density is corrected close to the symmetry axis and consequently the momentum and energy equations as well for a fluid system. We start with the same approach as proposed in [START_REF] García-Senz | Axisymmetric smoothed particle hydrodynamics with selfgravity[END_REF], where the corrected density equation is given by

a N b ab b a f W m 1 1       (29)
where a   is the new corrected 2D density for particle a and a f 1 is the correction function which is prominent close to the symmetry axis and equals 1 far away from the symmetry axis such that the density of the particle tends to zero as the radial distance of the particle goes to zero. A detailed derivation for the correction factor a f 1 can be found in [START_REF] García-Senz | Axisymmetric smoothed particle hydrodynamics with selfgravity[END_REF]. The function a f 1 for a particle a is given by Accepted Version of : Joshi, S., Franc, J.P., Ghigliotti, G. et However, the above corrective function is only applicable to the cubic-spline kernel function which is used in the present study. The above-mentioned density correction function is applicable to any axisymmetric SPH solver using a cubic-spline kernel. It applies to all materials irrespective of it being fluid and solid. However, the further derivations required for the momentum equations due to the density correction will vary depending on whether we solve for a solid or a fluid. The momentum correction derivation provided in [START_REF] García-Senz | Axisymmetric smoothed particle hydrodynamics with selfgravity[END_REF] is only applicable for a fluid solver and not a solid solver, we thereby present a detailed derivation of the corrected equation in the case of a solid solver.

                                       2 1 2 1 60
Since the density equation has been re-written as equation 29, the momentum equations should also be corrected to obtain correct acceleration and velocities for a particle close to the axis. For instance, the acceleration in the radial direction for a cylindrical system is given by

                r z r r dt dv rz rr rr r a       1 (31)
We start with the first term on the right hand side of equation 31:

r dt dv rr I r a      1 ( 32
)
This is calculated using the following identities:

                     rr rr rr rr rr rr r r r r r r r r       ) ( 1 ) ( (33
                                               2 2 1 ) ( 1 1 ) ( 1 (34)
Using corrected densities and the above identities, we can express equation 32 as

                  rr rr rr rr rr I r a r r r r r r dt dv 2 2 2 ) ( 2 2                            ( 35
)
where the derivative of the corrected density can be expressed as (36)

A similar approach is used for the third term on the right hand side of equation 32: which can be re-written as

r f r f r f r            ) ( ) ( )) ( . ( 1 1 1         Accepted
                               rz rz rz III r a r z z r z r dt dv 2 2 2 2 (38)
Finally, all the terms are added together to give the derivative of the radial velocity w.r.t. time

                        2 2 2 2 2 2 2                              rz rz rr rr r a r z z r r r r r dt dv (39)
The same approach is applied to derive a similar expression for derivative of z-velocity w.r. The above equation has the same form as equation 32 (the difference being that the quantities inside the derivative are different components of stress), hence applying the same identities as applied for 33 lead to

                                                    
             rz rz rz III z r r r r dt dv 2 2 2 2                (44)
Finally, all the terms are added together to give the derivative of the axial velocity (z-direction) w.r.t. time The resulting momentum equation given by equations 39 and 45 are rewritten in discrete SPH form as :

                                         
Accepted Version of : Joshi, S., Franc, J.P., Ghigliotti, G. et

f r r m dt dv 1 1 1 2 2 1 2 2 2 2 2 2                                                   (46)
1 2 2 1 2 2 2 2 2                                             (47)
The above momentum equations give a mathematically consistent solution for the treatment of particles close to the symmetry axis.

c. Validation test case for the axisymmetric SPH solid solver

The axisymmetric elasto-plastic solver is implemented and the simulation of an indentation is carried out as a test case to validate the solver. The indentation test is chosen for two primary reasons. Firstly, the indentation will induce a maximum displacement in vicinity of the symmetry axis and hence it will be a good test to identify the possible errors caused by the new scheme developed in the paper. Secondly, the load applied during indentation is chosen to be similar to that of cavitation due to a single bubble collapse. Figure 4 shows the domain used for indentation simulations. Stainless steel A-2205 is simulated, the material properties are as follows [START_REF] Roy | Determination of cavitation load spectra-Part 2: Dynamic finite element approach[END_REF]: Young's Modulus E=186x10 9 Pa, Poisson's ratio =0.3 and yield curve given by

)) log( 1 .( ) ( 29 . 0 0 p p Y Y K          (48)
where 0 Y  = 508x10 6 Pa, K = 832x10 6 Pa, α = 0.031. The ultimate plastic strain to failure is set to A localized indentation load is applied to the material using boundary particles with a prescribed velocity along the z-axis which depends on the distance to the axis. The velocity profile is shown in the right of Figure 4. Such a localized profile should create a continuously growing pit in the material. Consequently, it is expected that the material will first start to deform elastically, then some plasticity will develop in the most loaded region and finally some damage will initiate where the cumulated plastic strain exceeds the ultimate strain. The reason for using such a loading instead of an indenter itself is because using an indenter would require a free surface of the solid where indentation would be imposed. Free surfaces in solid or fluid SPH codes require special treatment as the particles near the free surface do not see a full kernel but a truncated one. To solve this most solvers, use density normalization on algorithms, these algorithm work quite well with 2D and 3D solvers where all the particle have the same mass. However, for an axisymmetric solver where mass of the particle is dependent on the initial radius and the mass in not constant for all particles, such a density normalization scheme does not work well specially near the axis where density calculations are not straightforward. Hence to avoid such issues a layer of boundary particle is used such that the top surface could also see a full kernel due to the presence of the boundary particles. A velocity type boundary condition is imposed to create a continuously growing pit.

To validate whether the scheme and the equations derived using the density correction (equation 29, 46, 47) are capable of simulating the solid behavior with an acceptable accuracy, the same case is setup in FEM and the axisymmetric SPH solver (with and without the novel mathematically consistent solution near the symmetry axis) is compared against FEM results. The FEM simulations have been performed with CAST3M [27] using 3456 8-node cells with a minimum mesh size of 2 μm and the same velocity profile is applied to the top nodes of the mesh as Dirichlet boundary conditions. The SPH simulation has been performed using 125570 particles and initial inter-particle distance of 2.5 μm and a smoothing length 2h = 6.5 μm. Firstly, we demonstrate that the existing method to solve solid mechanics in axisymmetric SPH without using the corrections derived in this paper (using equation [START_REF] Omang | SPH in spherical and cylindrical coordinates[END_REF]27,[START_REF] Pineau | An axisymmetric solid SPH solver with consistent treatment of particles close to the symmetry axis[END_REF] does not reproduce the solid behavior close to the axis with acceptable accuracy. Figure 5a shows comparison of stress components (left axis of the graph in figure 5a) and plastic strain (right axis of the graph in figure 5a) along line AB (dashed in red) shown in the left of Figure 4 when the pit depth has reached 6 μm. It is clearly observed that the stress and plastic strain near the axis are not well reproduced, a maximum error of 15% in stress components close to the axis and an error of 60% for plastic strain close to the axis is observed. To compare the two solvers, we plot all the stress components (left axis of the graph in figure 5b) and plastic strain (right axis of the graph in figure 5b) along line AB (dashed in red) shown in the left of Figure 4 when the pit depth has reached 6 μm. Figure 5a shows the comparison between SPH and FEM simulations. The average error is less than 1% suggesting an excellent agreement between the SPH and FEM results, hence validating the corrected density and momentum equations derived in this paper (equation 29, 46, 47). Also contours of plastic strain obtained from the FEM and SPH simulations are shown in figure 5c which shows that the FEM and SPH produce very similar plastic strain patterns in the material. a downward velocity with a Gaussian shape as plotted on the right (where vertical axis denote the Z-axis and the horizontal axis denote the radial direction). The extent of the velocity profile R is defined as the distance from the center where the velocity is 1% of the peak value, a non-uniform but constant velocity is given to the indenter.

SPH simulations for solid response a. Mass loss simulation and convergence study

We further use the solver to demonstrate that it is possible to run mass loss simulations using the axisymmetric SPH solver. We keep the same domain and the same simulation parameters as demonstrated in Figure 4. To account for mass loss, a damage criterion is to be introduced into the model. A thorough review of such models is provided by Pineau et. al. [START_REF] Pineau | An axisymmetric solid SPH solver with consistent treatment of particles close to the symmetry axis[END_REF],

Accepted Version of : Joshi, S., Franc, J.P., Ghigliotti, G. et these models are based on void growth theory wherein the void coalesce either under tension or shear to form cracks which eventually leads to damage (for details on these models please refer to the review article by Pineau et. al. [START_REF] Pineau | An axisymmetric solid SPH solver with consistent treatment of particles close to the symmetry axis[END_REF] and the references therein). However, in this paper, to demonstrate that it is possible to simulate mass loss simulation with SPH, we use a simple damage criterion based on the rupture strain [START_REF] Karimi | Phenomenological model for cavitation erosion rate computation[END_REF]. A critical value of rupture strain is defined to predict whether the SPH particle is damaged or not. If the plastic strain is above the rupture strain the material can sustain, u, the particle is assumed damaged. Therefore, once the plastic strain for a particle reaches a threshold value of rupture strain (εu = 0.03), the particle is considered as fully damaged and the particle thereafter does not contribute to the nearby particle interactions. This has been adopted similar to the approach in past studies where FEM cells are removed when a certain cell reaches a damage threshold [START_REF] Chandra | Modeling and analysis of material behavior during cavitation erosion[END_REF][START_REF] Roy | Determination of cavitation load spectra-Part 2: Dynamic finite element approach[END_REF]. However, since it is difficult to validate a mass loss simulation, it is not straightforward to verify whether the approach used for damaged in the present study is the best scheme. The value of rupture strain under tension is around 0.25 for stainless steel, however to demonstrate the capability of SPH to solve for damage and also keeping the simulation time within acceptable range, we chose a lower value of rupture strain as 0.03. Figure 6 shows a sequence of images from a simulation where such a damage criterion is used, the first image (Figure 6a) shows the plastic strain developing as the boundary particles start indenting the solid, the following image (Figure 6b) shows a higher plastic strain with peak plastic strain occurring at a certain depth below the top surface of the material. As expected the crack starts to propagate from the region of maximum plastic strain, the third image (Figure 6c) shows the crack propagating inside the material and finally (Figure 6d) material removal. However, it should be noted that mass loss simulations can be sensitive to the grid size in computational mechanics. Hence, we present a convergence study of the mass loss simulation while varying the interparticle distance. The simulation reported in figure 6 uses an interparticle distance of 2.5 µm, we conducted the same mass loss simulation for 3 more interparticle distances (2, 3.5 and 5 µm). In figure 7 we plot mass loss curves for all four inter-particle distances. The results show less than 4% difference between the highest and the lowest mass loss at the end of the simulations while the difference between 2 µm and 2.5 µm being less than 1%. Thus, it is verified that the mass loss simulation reported in the article are accurate enough and any increase in number of particle would add to the computational cost without any significant improvement in accuracy.

b. Parametric study for different extend of indent profile

A parametric study is performed to obtain mass loss curves for different extend of the applied velocity profile as indicated by radius R in Figure 4, keeping the velocity of indentation and other parameters the same. It should be noted that R is not an actual radius but just the extent of the velocity profile till it goes to 1% of the peak value. The mass loss curves vs. time are plotted in figure 11 for three different values of the extend of velocity profile (R) with coefficient C in equation 11 (Johnson-Cook model) varied as 0 and 0.031 to obtain non-strain rate sensitive response and strain rate sensitive response respectively. It can be noted that over a long period of time the larger radius of indent (R) is able to produce larger mass loss. This is expected, since a larger extent of the applied load is able to create plastic deformation in a larger volume of material, hence producing higher mass loss. However, what is intriguing is the lower incubation period for a smaller radius of indent (R). Incubation period is indicated as A in Figure 10. It is defined as the time required for material under load to initiate mass loss. The incubation period is an important parameter for most cavitation studies as it marks the initiation of material erosion. Cavitation erosion is roughly inversely proportional to the incubation time, which means that a longer incubation time would mean a slower rate of erosion, however the trend observed from the mass loss curves in figure 11 is quite the opposite. An important aspect to investigate would be to understand how does the smaller radius of indent (R) cause a quicker damage (low incubation time) in the material even though it produces higher rate of erosion. To investigate this, Figure 8 shows a comparison of the stress profiles and the plastic strain for two different R (0.2 mm and 0.4 mm) along the line AB shown in figure 4. It can be observed that close to the symmetry axis, all the stress components are higher for the R=0.2 compared to R=0.4 leading to a higher plastic strain and hence a smaller incubation period. However, the larder radius of indent (R=0.4) is able to plastify a larger volume of material and hence over a longer period of time shows higher mass loss. It should be noted that the type of loading simulated here (indentation generated by the indent velocity) is quite different from the impact load due to bubble collapse and hence it remains to be seen whether such a phenomenon (relatively higher incubation time for higher erosion rate) could be observed for a fluid-structure interaction simulation of a bubble collapse near a solid. If indeed a smaller bubble could lead to a lower incubation period and a low rate of erosion compared to larger bubble remains to be investigated in further studies. An important aspect of cavitation loading is that it can produce high strain rates and hence could lead to strain rate hardening in the material, thus lowering the rate of erosion. We numerically investigate the phenomenon by varying the coefficient C in equation 11 as 0 and 0.031, where 0 corresponds a non-strain rate sensitive material and 0.031 for a strain rate sensitive material like stainless steel. Rest of the material constants and simulation parameters are kept the same while C is varied from 0 to 0.031. We plot plastic strain for the two cases in figure 9 along the line AB (marked in figure 4), the strain rate insensitive simulation produces 10% higher plastic strain compared to the strain rate sensitive simulation. A similar behavior could be observed in figure 11 where mass loss shows a deviation of around 8-15% between a strain rate insensitive and sensitive simulation. Incubation time if found to be 3-8% higher for the strain rate sensitive case (C=0.031) compared to the strain rate insensitive case (C=0). However, it should be noted that the maximum strain rate during these simulations were observed to be of the order of 10 3 s -1 , numerical and experimental cavitation studies indicate that the strain rate involved in cavitation erosion process could be as high as 10 5 or 10 6 s -1 [START_REF] Roy | Determination of cavitation load spectra-Part 2: Dynamic finite element approach[END_REF][START_REF]Fundamentals of Cavitation[END_REF]. Hence the strain rate hardening effects in cavitation erosion can be expected to be much more prominent during cavitation loading. *Corresponding Author, Shrey Joshi: shrey.joshi@simap.grenoble-inp.fr Figure 11. Mass loss curves for different extend of velocity profile (R) on a stainless steel A-2205 specimen obtained from SPH simulations. The indentation depth is kept constant at 6 µm for all cases while the indenter radius is varied from 0.1 mm to 0.4 mm.

Conclusion and future work

An axisymmetric SPH solid solver is developed that can treat the particles close to the symmetry axis in a mathematically consistent way. The solver is capable of solving elasto-visco-plastic simulations with material damage and strain rate effects. An indentation test case is selected as a validation test case for reasons that firstly it can test the capability of the scheme in dealing with phenomenon near the symmetry axis and secondly that indentation loading is similar to cavitation loading. The results are compared against FEM results for the same case and an excellent agreement is observed. It is also demonstrated that the new formulation derived in the paper provides a much more accurate result near the symmetry axis compared to the existing method.

Mass loss curves were computed to demonstrate the capability of the model to solve for material damage. In this first study, the damage behavior was simplistic. A more precise and sophisticated damage model should depend on various aspects such as [START_REF] Pineau | An axisymmetric solid SPH solver with consistent treatment of particles close to the symmetry axis[END_REF]: cumulated plastic strain, stored energy, surface energy, stress triaxiality and strain rate. Various predictive models have been proposed so far in the past to determine fracture criteria considering all the above parameters. Future work will consist of including a more realistic damage models in the SPH code using experimental measurements.

The mass loss curves obtained for different extents of the applied load have shown that the smallest extent of the indenter velocity profile R is the fastest to initiate damage but produces the lowest rate or erosion The type of loading simulated here (indentation generated by the indent velocity) is quite different from the impact load due to bubble collapse. To understand the actual behavior, further this solver will be coupled to a fluid SPH solver to solve for cavitation bubble collapse near a solid, it remains to be seen if such a trend can also be observed when comparing erosion caused due to two different size of bubbles. Moreover, the present article only focuses on the cubic spline kernel and its correction factors, finding axiscorrections factors for higher-order interpolators could further enhance the axis-SPH simulations of indentation and cavitation problems.

Figure 1 .Figure 2 .

 12 Figure 1. Pressure-temperature phase diagram for water showing the two forms of vaporization in water, namely, boiling (state 1 to state 2) and cavitation (state 1 to state 3)

  relating the yield stress (σy) to the plastic deformation p  and the plastic strain rate p

Figure 3 (

 3 Figure 3(a) Sketch of the coordinate system and notation used to describe axisymmetric formulation, (b) shows a schematic of how particles are places in an axisymmetric system where real particles on right of the symmetry axis at a distance less than 2h from the axis are mirror imaged across to the axis as ghost particles.

  radial distance from the symmetry axis and a h is the smoothing length of the kernel.

  t. time (the first term of the right hand side of equation 26)

Accepted 8 *

 8 Version of : Joshi, S., Franc, J.P., Ghigliotti, G. et al. An axisymmetric solid SPH solver with consistent treatment of particles close to the symmetry axis. Comp. Part. Mech. (2020). https://doi.org/10.1007/s40571-019-00310-Corresponding Author, Shrey Joshi: shrey.joshi@simap.grenoble-inp.fr

Figure 4 .

 4 Figure 4. The figure on the left shows the computational domain for solid simulation, the boundary particles marked in blue are givena downward velocity with a Gaussian shape as plotted on the right (where vertical axis denote the Z-axis and the horizontal axis denote the radial direction). The extent of the velocity profile R is defined as the distance from the center where the velocity is 1% of the peak value, a non-uniform but constant velocity is given to the indenter.

Figure 5a . 8 *

 5a8 Figure 5a. FEM results compared against SPH (with ghost particles) without any correction near the axis (using equation[START_REF] Omang | SPH in spherical and cylindrical coordinates[END_REF] 27,[START_REF] Pineau | An axisymmetric solid SPH solver with consistent treatment of particles close to the symmetry axis[END_REF] for a pit of 6 microns in depth and 0.4 mm in radius for a stainless steel A-2205 specimen, the results are plotted on the horizontal red dotted line in figure4(0.12 mm below the top surface).

Figure 5c .

 5c Figure 5c. FEM and SPH (with correction scheme proposed in this paper) plastic strain contours for a pit of 6 microns in depth and 0.4 mm in radius for a stainless steel A-2205 specimen

Figure 6 .

 6 Figure 6. Sequence of images from SPH damage simulation of a stainless steel A-2205 specimen, (a) Shows the plastic strain accumulation in material without damage, (b) high plastic strain zone just beneath the top surface, (c) damage initiation from the point of highest plastic strain, (d) material loss due to a single indent (domain size shown in the figure is same as in figure 4).

Figure 7 .

 7 Figure 7. Mass loss vs time curves for simulations with four different interparticle distances: 2, 2.5, 3.5 and 5 µm.

Figure 8 .

 8 Figure 8. The graph shows components of stress and plastic strain plotted along the length of the red dotted line shown in figure 4 for two different values of extent of velocity profile R=0.2 (written as 0.2R in the graph) and R=0.4 mm (written as 0.4R in the graph). The left axis indicates values for all component of stress and the right axis indicates plastic strain along the length of the line.

Figure 9 .

 9 Figure 9. Shows plastic strain along the line AB marked for two cases, namely, one without strain rate sensitivity and the other without strain rate sensitivity.

Figure 10 .

 10 Figure 10. Characteristics of typical cumulative erosion versus exposure time curve. A = nominal incubation time; (B) = maximum erosion rate; (C) = terminal erosion rate; and D = terminal line intercept ([29]).

  al. An axisymmetric solid SPH solver with consistent treatment of particles close to the symmetry axis. Comp. Part. Mech. (2020). https://doi.org/10.1007/s40571-019-00310-8

	r dt dv	r	r	z	r	b	r

*Corresponding Author, Shrey Joshi: shrey.joshi@simap.grenoble-inp.fr

  al. An axisymmetric solid SPH solver with consistent treatment of particles close to the symmetry axis. Comp. Part. Mech. (2020). https://doi.org/10.1007/s40571-019-00310-8

*Corresponding Author, Shrey Joshi: shrey.joshi@simap.grenoble-inp.fr

  Version of : Joshi, S., Franc, J.P., Ghigliotti, G. et al. An axisymmetric solid SPH solver with consistent treatment of particles close to the symmetry axis. Comp. Part. Mech. (2020). https://doi.org/10.1007/s40571-019-00310-8

*Corresponding Author, Shrey Joshi: shrey.joshi@simap.grenoble-inp.fr

  al. An axisymmetric solid SPH solver with consistent treatment of particles close to the symmetry axis. Comp. Part. Mech. (2020). https://doi.org/10.1007/s40571-019-00310-8

*Corresponding Author, Shrey Joshi: shrey.joshi@simap.grenoble-inp.fr

  al. An axisymmetric solid SPH solver with consistent treatment of particles close to the symmetry axis. Comp. Part. Mech. (2020). https://doi.org/10.1007/s40571-019-00310-8

*Corresponding Author, Shrey Joshi: shrey.joshi@simap.grenoble-inp.fr

  Accepted Version of : Joshi, S., Franc, J.P.,Ghigliotti, G. et al. An axisymmetric solid SPH solver with consistent treatment of particles close to the symmetry axis. Comp. Part. Mech. (2020). https://doi.org/10.1007/s40571-019-00310-8

*Corresponding Author, Shrey Joshi: shrey.joshi@simap.grenoble-inp.fr

  Accepted Version of : Joshi, S., Franc, J.P.,Ghigliotti, G. et al. An axisymmetric solid SPH solver with consistent treatment of particles close to the symmetry axis. Comp. Part. Mech. (2020). https://doi.org/10.1007/s40571-019-00310-8

c. Strain rate effects

*Corresponding Author, Shrey Joshi: shrey.joshi@simap.grenoble-inp.fr

*Corresponding Author, Shrey Joshi: shrey.joshi@simap.grenoble-inp.fr

Acknowledgement

This project has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No 642536.

Authors also acknowledge the contribution and help from the Dr. Domingo García-Senz at Universitat Politècnica de Catalunya.

Appendix A

Kernel Function

The accuracy and stability of an SPH code is highly dependent on the choice of the kernel functions. They should satisfy several conditions such as positivity, compact support, and normalization. Also, ab W of a particle a should monotonically decrease with increasing distance from the particle and behave like a delta function as the smoothing length , tends to zero [START_REF] Monaghan | Smoothed particle hydrodynamics[END_REF][START_REF] Benz | Smoothed Particle Hydrodynamics: A review in The Numerical Modelling of Nonlinear Stellar Pulsations: Problems and Prospects[END_REF][START_REF] Liu | Mesh Free methods: Moving beyond the finite element method[END_REF]. The value of the Kernel function depends on the smoothing length h, and the non-dimensional distance between particles given by q = r/h, r being the distance between particles a and b. The parameter h, often called smoothing length, controls the size of the area/volume around particle a where contribution from the rest of the particles is considered, outside the radius h the contribution of the particles is either zero or negligible.

A cubic spline kernel has been used in the present work:

where D is 10/(7πh 2 ) in 2D and 1/(πh 3 ) in 3D.