Shrey Joshi

Jean Pierre Franc

Giovanni Ghigliotti

Marc Fivel

SPH modelling of a cavitation bubble collapse near an elasto-visco-plastic material

Keywords: Smoothed Particle Hydrodynamics, cavitation erosion, fluid structure Interaction, plastic strain, inertial effects, strain rate effects

A Smoothed Particle Hydrodynamics axisymmetric solver was developed in order to simulate the collapse of a single cavitation bubble close to an elastic-plastic material and study plasticity formation and hence material erosion. Findings indicate the relative importance of the material deformation due to the impact of the micro-jet and the shock wave that develop during collapse. A shock-wave dominated impact has a much higher material erosion ability compared to a micro-jet impact. Strain rate is found to have a significant effect on plastic deformation, with an overestimation of the plastic deformation up to 60% if strain rate effects are neglected in the case of stainless steel A2205. We also demonstrate that, although the impact pressure is maximum just below the collapsing bubble, maximum plastic strain occurs at a radial offset from the symmetry axis. This is the result of inertial effects that have an impact on both the magnitude and the position of the plastic domain in the material. A new non-dimensional parameter called effective pressure is introduced that can predict plastic strain location accurately for higher stand-off ratios. Alternatively, a characteristic time analysis also shows that it can be used for prediction of plastic strain zone in the solid for detached cavities.

Introduction

Understanding cavitation erosion better has been a constant quest for cavitation research. The simplest cavitation erosion empirical models relate the erosion resistance derived from a standard cavitation test [START_REF]Advanced Experimental and Numerical Techniques for Cavitation Erosion Prediction[END_REF][START_REF] Momma | A study of pressures and erosion produced by collapsing cavitation[END_REF][START_REF] Chahine | Noise and Erosion of Self-Resonating Cavitating Jets[END_REF][START_REF] Escaler | Cavitation erosiion prediction from inferred forces using material resistance data[END_REF][START_REF] Hammitt | Damage to solids caused by cavitation[END_REF][START_REF] Lee | Liquid impact erosion mechanism and theoretical impact stress analysis in TiN-coated steam turbine blade materials[END_REF] to the mechanical properties of the material such as Yield stress, Young's Modulus, hardness, toughness etc. Other kinds of empirical models relate the erosion resistance to the pitting rate or flow aggressiveness most often without considering material properties [START_REF] Hattori | Prediction method for cavitation erosion based on measurement of bubble collapse impact loads[END_REF][START_REF] Futakawa | Off-line tests on pitting damage in Mercury target[END_REF][START_REF] Soyama | Estimation of incubation time of cavitation erosion for various cavitating conditions[END_REF]. Models that can account for both the fluid and the solid behavior were developed [START_REF]Fundamentals of Cavitation[END_REF][START_REF] Karimi | Phenomenological model for cavitation erosion rate computation[END_REF][START_REF] Kato | Possibility of Quantitative Prediction of Cavitation Erosion Without Model Test[END_REF][START_REF] Berchiche | A Cavitation Erosion Model for Ductile Materials[END_REF], but some studies show that they fail to match the experimental results [START_REF] Karimi | Phenomenological model for cavitation erosion rate computation[END_REF]. Considering the complex phenomenon of cavitation and the various aspects not considered in these models, it is needed to have better ways to model material erosion. A lack of accurate numerical model for erosion calculation has been a bottleneck in understanding cavitation erosion. Most models as mentioned above either are empirical or phenomenological models that oversimplify the phenomenon of cavitation erosion. However, a few studies have tried solving for cavitation erosion using numerical methods such as FEM or even fluid structure interaction solvers that can provide a much better understanding. We focus on these studies and the shortcomings therein in order to improve the numerical scheme in the present work. Some studies have used an inverse approach based on the Finite Element Method where it is assumed that the impact pressure due to a single bubble collapse has a Gaussian distribution in space [START_REF] Roy | Determination of cavitation load spectra -Part 1: Static finite element approach[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 elastoplastic deformation of metallic materials[END_REF].

                - = 2 2 exp H H d r σ σ (1)
where, r is the radial distance to the symmetry axis and H σ is the maximum amplitude of the hydrodynamic impact pressure and H d is the diametric extent of the load. The inverse calculation provides the characteristics and the pressure distribution that could have led to a pit of given diameter and depth. This study was extended to dynamic loading [START_REF] Roy | Determination of cavitation load spectra-Part 2: Dynamic finite element approach[END_REF] by adding a temporal evolution to the pressure, modifying equation 1 to the following equation

                - -                 - = 2 max 2 exp 2 exp H H H t t t d r σ σ ( 2
)
Where t is time and tmax is the time when

H σ σ =
and tH is the characteristic impact rise duration. In a similar way dH represents the radial spread of the load. The authors found that at high frequency or small characteristic loading time, the pits formed for the same pressure and radial spread appear to be smaller than the pits for the static case i.e. without temporal evolution of pressure. The results were attributed to strain rate and inertial effects. However, it would be interesting to analyses these effects under the actual loading due to a bubble collapse which is much more complicated than a simple Gaussian pressure.

To understand cavitation erosion in entirety, some studies have focused on fluid structure interaction (FSI) solvers [START_REF] Hsiao | Modelling of material pitting from cavitation bubble collapse[END_REF][START_REF] Turangan | Numerical studies of cavitation erosion on an elastic-plastic material caused by shock-induced bubble collapse[END_REF]. An initial footstep towards a numerical understanding of cavitation using FSI would be to understand the cavitation erosion phenomenon for a single cavitation bubble. Hsiao et al. have attempted to develop a cavitation FSI solver wherein both fluid and solid response can be captured [START_REF] Hsiao | Modelling of material pitting from cavitation bubble collapse[END_REF]. In order to solve for the fluid bubble collapse, they use three codes: two fluid solvers, namely, boundary element method code (3DYNAFS-BEM) and a finite difference code (GEMINI), whereas the solid response was simulated using a finite element code (DYNA3D). The aim of having two fluid solvers is to create a hybrid scheme wherein the bubble collapse is solved using two different schemes which depend upon the kind of system to be solved. Cavitation bubbles create no shock waves during the growth phase and even during the collapse, shock waves are generated only towards the end of the collapse. Hence, they use the incompressible BEM until the end of the collapse where the flow is mostly incompressible. Once the shock wave generation is expected, the solver is switched to a compressible finite difference solver. At each iteration, the solid solver takes the pressure as input and solves for the material response. The code solves the material using a linear isotropic hardening law, however, most materials exhibit non-linear isotropic hardening. Also, the code does not consider strain rate hardening effects which are significant for fast dynamic phenomena like cavitation erosion.

Another work by Turangan et al. [START_REF] Turangan | Numerical studies of cavitation erosion on an elastic-plastic material caused by shock-induced bubble collapse[END_REF] studied fluid structure interaction using a Free-Lagrange scheme, FLM, that incorporates the compressibility, multi-phases and elastic plastic solid models. They simulated the collapse of 40 µmradius single bubbles attached to/near rigid and Aluminum walls by a 60 MPa-lithotripter shock, and the collapse of a 255 µm-radius bubble attached to a 25 µm-thick Aluminum foil by a 65 MPa-lithotripter shock. They simulated the bubble with different stand-off distance and found that an attached cavity can produce larger deformation compared to a detached cavity. However, again the material model does not consider non-linear hardening and strain rate sensitivity. A comprehensive cavitation FSI numerical model should ideally take into account the elasto-plastic behavior of the material as well as the development of damage. Since cavitation is a fast process, strain rate dependent hardening behavior should be included in the model. Also, a fluid structure interaction scheme is required. During bubble collapse, high intensity shock waves are produced along with micro jet. The shock wave travels through the fluid to the solid, a part of the wave is reflected back into the liquid and the rest gets transmitted to the solid [START_REF] Francois Axisa | Modelling of Mechanical Systems: Fluid-Structure Interaction[END_REF]. A two-way coupling is required to model this behavior properly.

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) technique, and 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. • Conventional methods of coupling a FVM solver to a FEM solver can be quite complicated and require 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 since both the fluid and the solid are solved using a Lagrangian formulation and hence such deformations are already taken care of via particle movement. 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]. The 2D solver is modified to 2D axisymmetric for the fluid solver. The axisymmetric solid solver is developed in-house. However, axisymmetric SPH solver suffers from an inconsistent definition of density near the symmetry axis. A new method to mathematically resolve this issue for the axisymmetric solid SPH solver inspired by an approach used in fluid axisymmetric SPH simulations is used in the present study. The axisymmetric fluid and solid solvers are validated against Rayleigh-Plesset collapse and FEM simulation respectively. The two solvers are eventually coupled together to obtain a fully coupled FSI solver capable of solving single bubble collapses over a solid medium to obtain elastic-plastic response using the Johnson-Cook model [START_REF] Johnson | Fracture characteristics of three metals subjected to various strains, strain rates, temperature and pressure[END_REF]. The paper presents simulations for a detached cavity and an attached cavity to understand the distinct material response the two cavities can produce. Attached and detached refers to whether the cavity surface is attached or detached to the solid medium respectively.

Methodology

A detailed methodology to change a 2D SPH solver to an axisymmetric SPH solver is provided in [START_REF] Joshi | An axisymmetric Solid SPH solver with consistent treatment of particles close to the symmetry axis: Application to cavitation erosion[END_REF]. SPH algorithm for solving a fluid and a solid together within a single code is provided in the following sections.

Fluid solver

A major issue with axisymmetric SPH is the treatment of particles that are close to the symmetry axis. Indeed, large errors in density and consequently other quantities can be observed with particles at distance less than 2h from the axis. The fluid axisymmetric SPH scheme implemented in the present study uses a corrective function for density calculation near the symmetry axis. Further, the density correction leads to a corrected momentum equation [START_REF] García-Senz | Axisymmetric smoothed particle hydrodynamics with self-gravity[END_REF].

The 2D density is calculated using the following equation,

a N b ab b a f W m 1 1 × = ∑ = η  (3)
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. Axisymmetric calculations use 2D density which is related to the actual 3D density by the following:

a a a r ρ π η 2 = (4)
Detailed derivation for the correction factor a f 1 can be found in [START_REF] García-Senz | Axisymmetric smoothed particle hydrodynamics with self-gravity[END_REF]. The function a f 1 for a particle a is given by Once the density is calculated, the pressure at each particle can be computed using an equation of state. In the present work, we use for liquid water the Tait equation of state that is given by the following,

         ≥ ≤ ≤       - + - + - ≤ ≤       + - + = - - - -
        -         = 1 0 1 γ ρ ρ B P (6)
where γ = 7 and

γ ρ / c B 0 2 0 1 =
, 0 ρ = 1000 kg m -3 is the reference density and 0 c is the sound speed at the reference density.

Since the density equation has been re-written as equation 3, the momentum equations should also be corrected to obtain correct acceleration and velocities for a particle close to the axis. The acceleration in the radial and axial directions for a cylindrical system is given by, Additional information regarding the equations 9 and 10 can be found in the work by Monaghan [START_REF] Monaghan | Artificial viscosity for particle methods[END_REF]. Once the particle acceleration is obtained from equations 7 and 8, the acceleration is time integrated to obtain the velocities of each particle and eventually particles are moved at each time step according to their velocities (for details see [START_REF] Gómez-Gesteira | User Guide for the SPHysics Code v2.0[END_REF]).

Π - × -         × + - = ∑ ∑ 1 1 1 2 2 2 2 2 η π η η π η π     (7) ( )

Solid solver

A detailed derivation of corrected density and momentum equations near the symmetry axis for a solid SPH axisymmetric solver is presented in [START_REF] Joshi | An axisymmetric Solid SPH solver with consistent treatment of particles close to the symmetry axis: Application to cavitation erosion[END_REF]. The 2D density is calculated using equation 3 which remains the same as for the fluid. Velocities from the last time step are used to calculate the strain rate which is given by,

1 1 ( ) ( ) 2 n ij i i j j ab ab a a b b a b a j i b a a dW dW m v v v v dx dx ρ ε =   = - + -     ∑  (11)
The strain rate can be time integrated to obtain the strain and eventually the stress using the following relation for an elastic solid,

kk ij ij ij µ ε λδ ε σ + = 2 (12)
Once the stress is obtained, the acceleration of each particle is calculated using the following

dW f r r m dt dv 1 1 1 2 2 1 2 2 2 2 2 2 × + -         × + +         × + = ∑ ∑ η σ π η σ π η σ η σ π η σ η σ π θθ       (13)
dW f r r m dt dv 1 1 1 2 2 1 2 2 2 2 2 × +         × + +         × + = ∑ ∑ η σ π η σ η σ π η σ η σ π      (14)
Detailed derivation of these equations can be found in [START_REF] Joshi | An axisymmetric Solid SPH solver with consistent treatment of particles close to the symmetry axis: Application to cavitation erosion[END_REF]. Further the accelerations can be used to obtain the new position of each particle in the same way as for the fluid. The above equations are only valid for an elastic solid. However plastic deformation is a key aspect in material response under cavitation load. To calculate plastic strain, the Johnson-Cook model is used (temperature effects are neglected in the paper),

)] ln( C ][ ) ( B A [ ) T , , ( * p n p p p y ε ε ε ε σ   0 0 0 1 + + = (15)
where, A0, B0, C0, n are material constants, the reference plastic strain-rate of the quasi-static test used to determine the yield and hardening parameters A0, B0 and n. Details regarding plasticity calculation and return mapping algorithm for non-linear isotropic hardening are provided in Appendix B.

Fluid structure interaction in SPH

Since the solid and the fluid solver are coded within a single code, all particles irrespective of them being solid or fluid particles are solved at the same time. However, a few modifications are required for the particles close to the interface. Consider the situation in figure 1, where particle a is a solid particle and its kernel contains fluid particles, one of them being particle b. While calculating density of the solid particle a, the contribution from the fluid particle (which has a lower mass) can lead to underestimation of the density in the solid and overestimation of density in the fluid. In other words, a numerical density diffusion would be observed near the interface unless corrective measures are taken. To achieve this and avoid any numerical density diffusion, the density equation is modified to the following,

0 0 1 1 b a a N b ab b a f W m ρ ρ η × × = ∑ =  (16)
Where 0 a ρ and 0 b ρ is the density of particle a and b respectively at the start of the simulation.

To transfer force across the interface while solving all the particles together, for any interaction of particle across the interface, the total force on a given particle a due to particle b can be resolved in the r-z direction, while considering there is no interface but a continuum across the interface.

j F i F F z r ab + =  (17)
However, since there exists a discontinuity at the interface and in the absence of viscous forces, the force parallel to the interface should be zero. Hence the force vector is multiplied by the vector normal to the interface to obtain the final force on the particle.

j n F i n F F z z r r Interface ab + =  (18)
The above equation along with the momentum equation is enough to satisfy the Neumann boundary condition to transfer the forces across the interface. As the particles across the interface are solved together within a single code and at the same time, the velocity at the interface always satisfies the Dirichlet boundary condition. However, to avoid any penetration and avoid disorder, XSPH [START_REF] Monaghan | On the problem of penetration in particle methods[END_REF] is used which takes into account not just the velocity of the particles itself but also the velocity of the nearby particles while moving it.

Validations

Fluid solver validation

To validate the fluid solver, a classic case of a spherical bubble collapse is carried out using SPH and compared against the analytical solution to the Rayleigh-Plesset equation. Figure 2 shows the domain used for a bubble collapse simulation in a large medium. A void bubble of radius 0.099 mm is placed in a spherical domain which is 7 times the bubble radius. A pressure wave of 90 MPa is generated by the wave maker whereas the initial pressure inside the bubble is 0 Pa. All the domain walls are given a non-reflective boundary condition to avoid any wave reflections hitting the bubble. The flow is treated as non-viscous and no surface tension forces are considered in the model. The wave maker shown in the figure is a set of repulsive particles which apply force on the nearby particles, hence acting as a pressure wave generator. Figure 3 shows good agreement between SPH results and the analytical solution of Rayleigh-Plesset equation, which validates the fluid solver. The slight increase in bubble radius at the end of the collapse is due to the high pressure shock wave generated at the center by the collapse. The shock wave travels outwards and this outward travelling shock leads to a post-shock velocity in the same direction as the shock. That leads to a weak rebound of the particle near the center. Also the bubble radius does not go to zero in the figure 3, that can be explained by the fact that since particles have to keep a certain distance from each other (the interpartcile distance), the minimum radius at the end of the collapse is dependent on the interpartcile distance and would tend towards zero as the number of particles used in the simulation are increased.

Solid solver validation

An indentation test is used to validate the SPH solid solver against an FEM solver. Figure 4 shows the domain used for indentation simulations. Stainless steel A-2205 is simulated; the material properties are given in Table 1. In the validation step, the strain rate sensitivity coefficient C0 will be taken as 0 so the behavior will be assumed to be strain rate insensitive. The same case is setup in FEM and the axisymmetric SPH solver, with the novel mathematically consistent solution near the symmetry axis. The FEM simulations have been performed with CAST3M [24] using 3456 8-node elements 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. Figure 5 shows a comparison of stress components and plastic strain along the line AB (marked in figure 3) for the SPH and FEM simulation. The comparison shows excellent agreement with average error less than 1% between the SPH and FEM. Figure 6 shows the plastic strain contours obtained for the SPH and FEM simulation.

Results and discussion

In this paper, we present two different cases of bubble collapse over a deformable solid medium, one with the bubble attached to the solid and the other with the bubble detached from the solid. The following sections would present the results for the two cases and we point out some interesting findings about the response of the material due to cavitation loading. The material parameters correspond to A2205 duplex stainless steel including strain rate sensitivity effects as identified by Roy et al. [START_REF] Roy | Determination of cavitation load spectra-Part 2: Dynamic finite element approach[END_REF], the material properties are listed in Table 1. Figure 7a shows the sketch of the initial domain used for simulating a detached cavity collapse over a solid. A bubble radius of 0.15 mm and a standoff ratio (SR) of 1.33 are used for the initial domain (where SR is defined as D/R, D is the distance from the bubble center to the interface and R is the radius of the bubble). The figures 7b-d show a sequence of images for the collapse of the cavity and subsequent response of the material. To initiate the collapse, a pressure wave of 50 MPa is generated using a wave generator from the top of the fluid domain. Figure 7b shows contours of axial velocity (z-direction) in the fluid and plastic strain in the solid at τ = 0.92 (where τ is the simulation time at that instance divided by time required for the bubble to collapse), with the formation of a micro jet after the pressure wave hits the bubble (the dotted circle in figure 7b shows the initial bubble before the start of collapse).

Single bubble collapse: Detached cavity vs Attached cavity

Figure 7c shows contours of pressure in the fluid and plastic strain in the solid at τ = 1.053. A shock wave is generated due to the micro jet hitting the liquid on the other side of the bubble as shown in figure 7b. Figure 7c shows the shock just reaching the interface. This shock wave produces plastic deformation in the material as shown in Figure 7d. A maximum plastic strain of 0.0023 can be observed just beneath the interface and not at the solid surface. Surprisingly, the maximum plastic strain does not occur right beneath the center of the bubble but at an offset from the symmetry axis marked as distance dPM in figure 7d (dPM is around 1.1 times the initial bubble radius for the case presented in figure 7). This material behavior for the detached cavity will be dealt with in detail in section 4.4. 8b hits the material and produces a very tiny zone of plastic strain (max plastic strain 0.0042) at the center (the contours in figure 8c are magnified along the white dotted rectangle shown in figure 8b), it also demonstrates the formation of a toroid bubble after the micro jet impact on the solid, (d) plastic strain contours in the solid at the end of the simulation, contour in the material shown in figure 8d is magnified along the white dotted line in figure 8c, two zones of plastic strain can be seen: one due to the torioid collapse and the other at the center due to the micro jet and also the shock wave from the toroid collapse.

Figure 8a shows the sketch of the initial domain used for simulating an attached cavity collapse over a solid. A bubble radius of 0.15 mm and a standoff ratio (SR) of 0.66 are used for the initial domain. To initiate the collapse, a pressure wave of same amplitude 50 MPa is generated, still using a wave generator from the top of the fluid domain. As for the previous case of a detached cavity, a micro jet can be observed after the pressure wave hits the bubble (Fig. 6.8b).

Figure 8c shows contours of axial velocity (z-direction) in the fluid and plastic strain in the solid (the contours are magnified along the white dotted rectangle marked in figure 8b). A tiny zone of localized plastic deformation (maximum plastic strain 0.00042) at the center is observed due to the micro jet hitting the material. Also, a toroidal bubble can be observed towards the right of the micro jet. The toroidal cavity then collapses to further produce plastic deformation in the material as shown in Figure 8d corresponding to the end of the simulation, after unloading. As the toroid collapses, it produces plastic deformation below it (marked as plastic strain due to toroid collapse in figure 8d). Also, the shock wave due to the collapse of the toroidal bubble travels towards the center of the material and focuses on the axis to produce additional plastic deformation. A maximum plastic strain of 0.0042 is produced due to the micro jet and thereafter an additional plastic deformation is observed at the center due to the shock wave produced by the toroid collapse to give a final maximum plastic strain of 0.0083 at the center. This can be referred to as the primary zone of plastic deformation. A secondary zone of plastic deformation is observed below the point of toroid collapse with maximum plastic strain of 0.002.

Comparing the two cases, the results suggest that the micro jet generated by the collapse of an attached cavity (see fig [START_REF] Lee | Liquid impact erosion mechanism and theoretical impact stress analysis in TiN-coated steam turbine blade materials[END_REF].8c, before the toroid collapse) has an ability to cause a larger (almost 2 times larger) plastic strain in the material as compared to the shock wave generated by the collapse of a detached cavity, for the same magnitude of pressure wave initiating the collapse and the same bubble size. As a result, it could be conjectured that, the repeated collapse of an attached cavity would lead to a smaller incubation time compared to the repeated collapse of a detached cavity. But it should also be noted that the volume of material that is plastically deformed in case of a micro jet is miniscule compared to a shock wave impact (almost 800 times smaller). This would imply that even though the incubation time for material erosion might be lower for a micro jet collapse, the shock wave can plastify a much larger volume of material and hence the erosion rate should be higher for a shock wave impact. Hence it could be inferred that the material erosion ability of a shock wave is much higher than that of a micro jet.

Detached cavity material response

We get our attention back to the unexpected material response observed for the detached cavity. As shown in Figure 7d, the maximum plastic strain does not occur right below the center of the bubble but at an offset dPM from the symmetry axis. To understand this behavior, we plot the pressure in the fluid along the interface (r-axis as shown in figure 7a) in figure 9 at different times. The plot shows that the maximum pressure decreases as we move ahead in time (time 1 to time 6), which is the expected behavior for a spherical wave. Although the maximum pressure for all times is observed at the symmetry axis, the maximum plastic strain occurs at an offset dPM from the symmetry axis as shown by the red circle markers in figure 9. This behavior is contrary to the popular belief that cavitation erosion occurs where the maximum pressure is observed at the interface. Two hypotheses are considered to explain this behavior. The first hypothesis refers to a material effect by strain rate hardening, the second refers to dynamic effects in the relation between loading and plastic strain. These hypotheses are worked out in the following section.

respectively, where τ is the simulation time divided by the bubble collapse time), the right axis corresponding to the red circle markers refers to the plastic strain plotted along a line located just below the interface and parallel to it ( the line passes through the point of maximum plastic strain in the material).

Detached cavity: Strain rate effect

It could be due to strain rate effects, as the material can demonstrate hardening when strain rates are high. In order to examine this hypothesis, we eliminate strain rate dependency of yield stress by substituting C0 = 0 in equation 15, instead of C0 = 0.031. The contours of plastic strain for C0 = 0 are shown in figure 10. A maximum plastic strain of 0.0037 can be observed which is 60% higher than for the C0 = 0.0031 case. However, the offset dPM still remains the same for the two cases. Thus such a material behavior is not due to the strain rate effects.

Dynamic loading in cavitation

Since strain rate effects cannot account for the unexpected material behavior for the detached cavity case, it could well be the fast dynamics during cavitation loading that can lead to such a behavior. In order to analyze whether this is true, we follow the same strategy (as for the analysis of strain rate effects) by eliminating dynamic effects from the simulation. We achieve this by using the same load as seen by the material during an FSI simulation, but run a static simulation using an FEM solver since the present SPH solid solver cannot be used for static simulations.

Static vs dynamic loading in cavitation

Figure 11. Contour of plastic strain in the material for static simulation, simulated using FEM with pressure as boundary condition (pressure obtained from the interface pressure in SPH simulation at 500 time steps).

In order to setup a static case in FEM using the same load as experienced by the material in the FSI simulation, we use the load data between the time when the shock wave hits the solid material as shown in figure 7c (where the shock wave just reached the interface at its left end, let's call this instance T1) to when the shock wave travels along the interface and reaches the right end of the interface, let's call this instance T2. Between these times T1 and T2, we capture the pressure along the interface in the fluid from SPH simulation reported in figure 7, for 500 intermediary steps. Thereafter, the FEM solver [24] is used to setup a solid simulation with, as a boundary condition, the pressure data obtained at 500 time steps from the FSI SPH simulation. The minimum mesh size is set as 5 μm (using a QUA8 mesh type) and 5920 cells are used for the FEM simulation. The material properties are set the same as in the SPH simulation with α = 0 which refers to no strain rate effect. The FEM solver solves the material for each of the 500 steps one after the other. The plastic strain is accumulated in the material for all the 500 steps to obtain a static response of the material for the same load as the solid experiences during the SPH simulation. Figure 11 shows the response of the solid when subjected to a static loading simulated using FEM. Contours of plastic strain show that the static response does not have an offset dPM for the maximum plastic strain which occurs at the symmetry axis. This result would suggest that the reason for the offset of maximum plastic strain as shown in figure 7d is the dynamic loading. A dynamic simulation differs from a static one in terms of the rate of loading and unloading. A very fast rate of loading and unloading (dynamic simulation) would not allow sufficient time for the material to respond to the load because of inertial effects which are dealt with in the next section.

Inertial effects on material response

To understand the inertial effects on material response, we now look at the dynamics of the shock wave produced by the bubble collapse when it interacts with the fluid-solid interface. Figure 12 shows the shock wave interaction at the interface. The contours are taken from the SPH FSI simulation for detached cavity with SR = 1.33 as previously presented in figure 7. Contours of pressure are plotted in the fluid and plastic strain in the solid at a certain time during the simulation after the shock wave has hit the interface. A primary shock wave can be observed in the fluid which is generated by the bubble collapse. Once this shock wave reaches the interface, a part of the shock is reflected back in the fluid which is marked as a reflected shock wave in figure 12. The point of intersection of the primary and reflected shock wave corresponds to the maximum pressure along the interface as seen in the contours in figure 12.

Figure 12. Contours of pressure in the fluid and plastic strain in the solid, shows the primary shock wave (generated by the bubble collapse) and a reflected shock wave (generated by the reflection of primary shock wave at the interface).

In order to understand the phenomenon better, a schematic is presented in figure 13, demonstrating the shock wave evolution wrt time. Consider a spherical shock wave generated in the fluid at T=0 as marked in the figure 13 (T corresponds to a fictitious time used to exhibit the time evolution of the shock wave). The red dot corresponds to the point where the shock wave initiates, supposed to be at an axial distance αz from the interface. The shock wave then travels spherically outwards from the point of shock generation. At time T=3, the shock just reaches the interface. At T=4, along with the primary shock wave, a reflected shock wave is observed as the primary shock wave gets reflected from the interface. The point of intersection of the primary and reflected shock waves is marked as A, B & C corresponding to the time T=4, T=5 & T=6 respectively. These points correspond to the location of the maximum pressure along the interface at a given time as shown in figure 9. We will refer to this point where the primary and reflected shock wave intersect as the 'shock front' from now on in the article. Consider an SPH particle right below point A. This solid particle would respond to the pressure it experiences from the fluid and deform accordingly. If the shock front moves really quickly, the material might not have enough time to react to the load and deform as much as it would have had for the same magnitude of pressure under a static load. It is then of primary importance to compute the speed of the shock wave along the interface.

Figure 13. Schematic of time evolution of the primary shock wave and the reflected shock wave (times T=0 to T=6 correspond to fictitious time separated by a constant time interval, used to exhibit the position of the wave wrt time).

The location of the primary shock wave on the interface (as shown in schematic in figure 13) in an r-z plane can be represented by a circle and is given by, ( )

2 2 2 ) (t R z r z = - + α (19) where ) (t R
is the radius of the primary shock wave as a function of time. It can be written as

∫ = dt t u t R ) ( ) ( , ) (t u
is the radial velocity of the primary shock wave as a function of time, r is the radial position and z is the axial position of the wave front. Substituting z = 0 in equation 19 would give the position of the shock front (for example point A, B or C in figure 13) wrt time along the interface, derivative of the radial position gives the velocity of the shock wave along the interface (𝑣𝑣 𝑆𝑆𝑆𝑆 ) given by equation 20.

𝑣𝑣 𝑆𝑆𝑆𝑆 = 𝑐𝑐 𝑙𝑙 � 1 + � 𝛼𝛼 𝑧𝑧 𝑟𝑟 � 2 (20)
To obtain an analytical solution, we substitute 𝑐𝑐 𝑙𝑙 = 1500 m/s (speed of sound in water) and z α = 0.09 mm, which corresponds to the position where the micro jet hits the bubble surface for the detached cavity case with SR = 1.33. The solution of equation 19 & the equation 20) gives the position and velocity of the shock front wrt time. Figure 14 shows the plot of the theoretical velocity of the shock front wrt to the radial distance from the axis (legend marked as SW Intersection Velocity Analytical). The velocity of the shock front is high close to the symmetry axis and goes down till it reaches the speed of the shock wave itself, which is here set as 1500 m/s. The high velocity near the symmetry axis is due to the spherical shock front hitting a flat surface. As a result, the rate of loading and unloading is also high close to the symmetry axis, which does not allow the material sufficient time to respond to the load, hence exhibiting lower plastic strain near the symmetry axis even though the pressure acting on the material is the highest at the symmetry axis.

To sum up the argument, we have plotted the pressure and the velocity corresponding to the shock front wrt to the radial distance from the axis in figure 14 as obtained from SPH simulations. The computed velocity of the shock front obtained from SPH compares well with the analytical solution. The difference in values could be due to the fact that the shock wave velocity is considered as a constant (1500 m/s) for the analytical solution, whereas it actually depends on the pressure difference between the shock front and rear. The background contour (plastic strain) on the graph is set to scale with the x-axis of the graph which corresponds to the distance along the interface. It can be clearly seen that the pressure at the shock front (which is also the maximum pressure for a given instance) is maximum at the symmetry axis and decreases as we move along the interface. Moreover, the point at which the pressure decreases below the material yield stress (508 MPa) also corresponds well to the end of the plastic zone. Clearly, the background contour demonstrates that the maximum plastic strain occurs not at the point of maximum pressure but at an offset somewhere between 0.1 mm to 0.2 mm from the symmetry axis along the interface.

Figure 14. Plot for the velocity of the shock front (right axis) wrt its position along the interface obtained from both the SPH simulation and analytical expressions (equation 19). Plot also shows the pressure of the point of intersection (left axis) wrt its position along the interface obtained from SPH simulation The background contours

are for plastic strain in the solid and are set to scale with the x-axis.

Effective Pressure

A physical explanation was given in the last section that demonstrates why the plastic deformation occurs at a radial offset. However, it would be interesting to understand the phenomenon from a quantitative perspective. Moreover, in this section we aim at providing a criterion that could predict the zone where plasticity will occur based on quantities available in classic CFD calculations, namely pressure and velocity in the fluid. We define a new variable named effective pressure as given below,

( ) 2 0 l SW Y effective c v ) P ( P - - = ρ σ ( 21
)
where P is the pressure at the interface where the effective pressure is being calculated, 0 Y σ is the initial Yield stress for the material (508 MPa in the present simulation), ρ is the density of the solid , SW v is the velocity of the shock front along the interface at the point where the effective pressure is calculated and l c is the speed of sound in the liquid. The effective pressure is a non-dimensional parameter that should indicate the pressure that effectively leads to plastic strain in the solid. The numerator gives the excess pressure above the Yield stress, as only the load above the Yield stress could lead to plasticity. The denominator takes into account the inertial effects due to the density of the solid and the excess velocity of the shock front above the sound speed of liquid. The excess velocity is the shock front velocity ( SW v ) minus the speed of sound in liquid ( l c ). The term ( )

l SW c v -
makes sure that if the velocity of shock front along the interface is much higher than the sound speed, the effective pressure acting on the material (due to inertial effects) is reduced via the term ( )

2 l SW c v - ρ
. The density term in the denominator plays an important role since the phenomenon is linked to the inertial effects in the material, however for the present simulation the density effects can be neglected as we simulate the same material for different cases. .2 shows a maximum plastic strain at an offset which is due to the inertial effects as explained earlier. However, close to the symmetry axis a tiny zone of relatively high plastic deformation can also be observed, which is actually due to the impact of the micro jet. The effective pressure definition (equation 20) being based on the assumption that all energy that reaches the interface is via the shock wave impact, the effective pressure near the symmetry axis for SR=1.2 does not correlate with the micro jet impact. The lower the stand-off ratio, the higher the dynamic pressure experienced by the interface would be.

Hence the effective pressure definition can be used to predict plastic strain when the stand-off ratio is sufficiently large (SR >1.3). For lower stand-off ratios, the plastic strain close to the symmetry axis is under predicted by effective pressure and an alternative approach should be proposed that is the aim of the next section. An alternative approach to quantitatively analyze the material behavior is presented in this section. We focus on two quantities: material characteristic time (or characteristic response time) and loading characteristic time. We use the definition of material characteristic time as mentioned in Roy et al. [69], which is given by the depth of the plastic zone divided by the speed of sound in solid and which signifies the amount of time required for the material to react to a load. Figure 17a shows the depth of the plastic zone at various points along the interface from which the material characteristic time was computed. On the other hand, the loading characteristic time for a point at the interface is defined as the time during which the pressure at that point first rises above 100 MPa (black dashed horizontal line in figure 17b) and finally drops down below 508 MPa (initial yield stress for the material marked as green dashed line in figure 17b).The exceedance limit of 100 MPa is chosen in our case since a sharp change in slope of the pressuretime curve is observed around 100 MPa as evident from figure 17b, it is thus assumed that the loading due to bubble collapse starts at the time when the pressure reaches 100 MPa. It should be noted that if the material characteristic time is greater than the loading characteristic time, the material behavior would be dominated by the inertial effects. Figure 18 shows the plot of both material characteristic time and loading characteristic time along the interface for SR = 1.33, the contours of plastic strain are set as the background in the plot. It can be observed that for the region close to the symmetry axis where the material exhibits relatively low plastic deformation (region marked as A in the figure), the loading characteristic time is lower than the material characteristic time and for region B where the material characteristic time is lower than the loading characteristic time, the material has time to deform under the load and hence exhibits relatively higher plastic deformation

For the case of SR = 1.2, where the effective pressure fails to capture the plastic zone close to the symmetry axis, it would be interesting to see whether characteristic time plot can capture it better. Figure 19 shows the characteristic time plots for SR = 1.2. It can be observed that in region A where the plastic strain is relatively high, the material characteristic time is lower than the loading characteristic time thus predicting well the effect of inertia and dynamic pressure. As we move further along the interface into region B, the material characteristic time becomes marginally larger than the loading characteristic time, thus leading to a slightly lower plastic strain zone as the inertial effects come into play. Further, in region C the material characteristic time is much lower than the loading characteristic time thus the material exhibits higher strain rate since the inertial effects are not dominant. Thus, characteristic time analysis can be used to predict plastic strain zone in the material. However, the definition of the material characteristic time used in this study requires knowing the plastic strain depth and hence cannot be estimated from an only fluid simulation. Background image is the contours of plastic strain for SR = 1.2.

Experimental comparison

A well-known experimental study by Phillip and Lauterborn [START_REF] Philipp | Cavitation erosion by single laser-produced bubbles[END_REF] report findings for bubble collapsing over soft aluminum and the consequent damage that occurs due to one or more collapses. One of the results in their paper show damage patterns for standoff ratio 1.4 (which is close to the standoff ratio of 1.33 for detached cavities used in this study), where a ring shaped pattern is observed. The ratio of the damage ring radius to the maximum bubble radius is observed to be around 1.3 by Phillip and Lauterborn [START_REF] Philipp | Cavitation erosion by single laser-produced bubbles[END_REF] in the experimental findings which is close to the ratio of dPM and bubble radius of 1.1 for the numerical simulations in our case. Phillip and Lauterborn [START_REF] Philipp | Cavitation erosion by single laser-produced bubbles[END_REF] associate the occurrence of ring type damage to the fact the rebound bubble collapses over the damaged ring. However, it cannot be denied that it can be first the plastic strain developed due to the shock wave from the primary collapse (as shown by our simulation) and then additional plastic strain due to the micro bubbles in the rebound bubble that causes the damage. Thus even though there does not exist a direct evidence of the fact that shock wave causes a ring shaped deformation, the experimental and the numerical results suggest that it is very much possible that the ring can be due to both: first the shock wave and then the rebound bubble, this of course remains a topic of further investigation.

Conclusions and future work

An axisymmetric SPH FSI solver has been developed, which is capable of computing the collapse of a single bubble over an elastic-plastic solid following Johnson-Cook yield criterion. The fluid solver was validated against the analytical solution of the Rayleigh-Plesset equation and the solid solver against an FEM solver for an indentation case.

Simulations of the collapse of a detached and an attached cavity suggest that the micro jet generated during the collapse of an attached cavity has an ability to cause a larger maximum plastic strain in the material as compared to the shock wave generated during the collapse of a detached cavity. It is observed that for the same magnitude of pressure wave initiating the collapse and the same size of the bubble, the micro jet can produce twice the maximum plastic deformation compared to a shock wave. Hence, in case of repeated collapse, a micro jet dominated impact would exhibit a smaller incubation time compared to the detached cavity. On the other hand, the volume of material that is plastically deformed in case of a micro jet is miniscule compared to a shock wave impact (almost 800 times smaller). This would imply that, even though the incubation time for material erosion might be lower for a micro jet collapse, the shock wave can plastify a much larger volume of material and hence the erosion rate should be higher for a shock wave impact. Hence it could be inferred that the material erosion ability of a shock wave dominated collapse is much higher than that of a micro jet dominated collapse.

An important aspect of cavitation loading is the high strain rate associated to the phenomenon. The present results show that the strain rate effects can significantly affect plastic deformation in the solid since non-strain rate sensitive simulations produce around 61% higher plastic deformation for a detached cavity compared to a strain rate sensitive model in the case of stainless steel A2205.

An important and novel finding in the present article is the response of the material for a detached cavity. It is shown that the maximum plastic deformation does not occur at the center of collapse but at an offset from the center (cf. figure 7d). Even though the pressure experienced by the material is the highest at the center, it does not produce the maximum plastic deformation there. This phenomenon is due to inertial effects, since the material tends not to respond to the load as the rate of loading and unloading is extremely high. The effect is linked to the high velocity of shock front along the interface close to the symmetry axis. The study clearly demonstrates that maximum pressure does not always correspond to the location of maximum plastic deformation.

A new parameter called effective pressure is defined in the article which matches well with the plastic strain in the material. for a standoff ratio SR=1.33. Such a parameter can easily be used in CFD calculation to predict the location where plasticity will be localized. It is found that the effective pressure can also predict the location of the plastic zone for a higher standoff ratio SR=1.5. However, the definition of effective pressure is based on the shock wave impacting the solid and hence is only valid for large stand-off ratios where plasticity is governed by wave propagation only. For lower stand-off ratios the effective pressure under predicts the plastic strain near the symmetry axis which is due to the impact of the micro jet and not due to the shock wave.

Another approach to quantifying material response while considering inertial effects is then proposed based on a characteristic time analysis. The material characteristic time and the loading characteristic time are compared for two cases (SR = 1.2 & 1.33), where a lower material characteristic time relative to the loading characteristic time signifies more dominant inertial effects and hence relatively lower plastic strain. It is found that the characteristic time analysis predicts the zones of plastic deformation quite well. However, the definition of material characteristic time is such that it cannot be known without solving for the solid. And the equivalent plastic strain is defined as, Where t VM σ is the trial von-Mises stress in the material, P ε ∆ is the incremental plastic strain, 1 -N P ε is the plastic strain at time step N-1 and Y σ is the yield stress in the material. The above equation can be solved using numerical methods such as Newton-Raphson. Yield stress can be calculated using equation 48, however, the yield stress not just depends on the plastic strain but also on plastic strain rate. To calculate plastic strain, rate the following equation is used:

Acknowledgement

r

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

  -dimensional effective plastic strain rate) where 0 p ε is

Figure 1 .

 1 Figure 1. Interaction of solid and fluid particles near the interface, red circle represents the kernel of particle a.

Figure 2 .

 2 Figure 2. Simulation domain for bubble collapse simulation in 2D axisymmetric.

Figure 3 .

 3 Figure 3. Comparison of bubble radius vs. time, SPH against Rayleigh-Plesset solution for a bubble collapse simulation (domain shown in figure 2).

Figure 4 .

 4 Figure[START_REF] Escaler | Cavitation erosiion prediction from inferred forces using material resistance data[END_REF]. The figure on the left shows the computational domain for solid simulation, the boundary particles marked in blue are given a downward velocity with a Gaussian profile with the distance to the axis shape as plotted on the right. 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 5 .

 5 Figure 5. FEM results compared against SPH (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 figure 3 (0.12 mm below the top surface).

Figure 6 .

 6 Figure 6. FEM and SPH 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 7 (

 7 Figure 7(a) Sketch of the initial simulation domain for axisyymmetric SPH FSI simulation of a detached cavity of radius R = 0.15 mm and standoff ratio SR = 1.33, (b) shows contours of axial velcosity in the fluid and plastic strain in the solid, a micro jet formation can be observed (the dotted semi-circle represents the initial bubble shape)at τ = 0.92, (c) contours of pressure in the fluid and plastic strain in the solid, a shock wave generated due to the bubble collapse can be observed as it just recahes the interface at τ = 1.053, (d) plastic strain contours in the solid at the end of the simulation, the maximum plastic strain occurs at a radial offset dPM from the symmetry axis.

Figure 8 (

 8 Figure 8(a) Sketch of the initial simulation domain for axisyymmetric SPH FSI simulation of an attached cavity of radius R = 0.15 mm, (b) shows contours of axial velocity in the fluid and plastic strain in the solid, a micro jet formation can be observed, (c) contours of axial velocity in the fluid and plastic strain in the solid, the micro jet shown in figure8bhits the material and produces a very tiny zone of plastic strain (max plastic strain 0.0042) at the center (the contours in figure8care magnified along the white dotted rectangle shown in figure8b), it also demonstrates the formation of a toroid bubble after the micro jet impact on the solid, (d) plastic strain contours in the solid at the end of the simulation, contour in the material shown in figure8dis magnified along the white dotted line in figure8c, two zones of plastic strain can be seen: one due to the torioid collapse and the other at the center due to the micro jet and also the shock wave from the toroid collapse.

Figure 9 .

 9 Figure 9. Plot of pressure (left axis) vs the distance along the interface for different times: Time 1 to Time 6 ( Time 1-6 refers to τ =1.053, 1.106, 1.16, 1.2, 1.24 & 1.28 respectively, where τ is the simulation time divided by the bubble collapse time), the right axis corresponding to the red circle markers refers to the plastic strain plotted along a line located just below the interface and parallel to it ( the line passes through the point of maximum plastic strain in the material).

Figure 10 .

 10 Figure 10. Contour of plastic strain in the material for C0 = 0 in equation 15 (refers to the case with no strain rate effect).

Figure 15 .

 15 Figure 15. Contours of plastic strain for 50 MPa driving pressure and 0.15 mm bubble radius for (a) Stand-off ratio SR = 1.2, (b) stand-off ratio SR = 1.33 & (c) stand-off ratio SR = 1.5. We present three cases with different stand-off ratios (1.2, 1.33 & 1.5). The contours of plastic strain are shown in fig 6.15 a-c. It can be observed that the three cases produce distinctly different patterns of plastic strain. The effective pressure given by equation 20 will be used for all the three cases to check whether the effective pressure can qualitatively match the plastic strain in the solid, thus giving us a parameter (effective pressure) that can be used in fluid simulations to indicate the erosion prone areas better while considering the inertial effects.

Figure 16 .

 16 Figure 16. Plot of plastic strain and effective pressure for (a) Stand-off ratio 1.2, (b) stand-off ratio 1.33 & (c) stand-off ratio 1.5 for 50 MPa driving pressure and 0.15 mm bubble radius. Plastic strain in the material (plotted along a line parallel to the interface in the material passing through the maximum plastic strain position) and effective pressure at the interface are plotted for all three cases (shown in fig 6.15) in fig 6.16 a-c. It can be observed that effective pressure defined by equation 20 qualitatively matches well with plastic strain in the material for SR = 1.33 & 1.5. However, for SR = 1.2, the plastic strain does not match equally well with effective pressure close to the symmetry axis. The contour in fig 6.15(a) for SR =1.2 shows a maximum plastic strain at an offset which is due to the inertial effects as explained earlier. However, close to the symmetry axis a tiny zone of relatively high plastic deformation can also be observed, which is actually due to the impact of the micro jet. The effective pressure definition (equation 20) being based on the assumption that all energy that reaches the interface is via the shock wave impact, the effective pressure near the symmetry axis for SR=1.2 does not correlate with the micro jet impact. The lower the stand-off ratio, the higher the dynamic pressure experienced by the interface would be.Hence the effective pressure definition can be used to predict plastic strain when the stand-off ratio is sufficiently large (SR >1.3). For lower stand-off ratios, the plastic strain close to the symmetry axis is under predicted by effective pressure and an alternative approach should be proposed that is the aim of the next section.

4. 4 . 4 .

 44 Figure 17. (a)Contour of plastic strain showing depth of plastic strain zone at different point on the interface. The material characteristic time is defined as the depth of the plastic strain divided by the speed of sound in the solid, e.g. material characteristic time at point P1 can be calculated as d1/cs (where cs is the sound speed in the solid), (b)shows a typical pressure vs time curve for a certain point on the interface, to estimate the loading characteristic time we consider the time from when the pressure is above 100 MPa till it drops down to 508 MPa (initial yield stress).

Figure 18 .

 18 Figure 18. Plot shows material characteristic time and loading characteristic time along the interface for SR= 1.33, plot background image is the contours of plastic strain for SR = 1.33.

Figure 19 .

 19 Figure 19. Plot shows material characteristic time and loading characteristic time along the interface for SR= 1.2.Background image is the contours of plastic strain for SR = 1.2.

2D density of particle a a η 

 η 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. Nomenclature W SPH Kernel function ha kernel smoothing length Pa pressure of particle a ma mass of particle a ρa density of particle a ρa0 Initial density of particle a (at the start of the simulation) σ stress tensor ε strain tensor i a v velocity of particle a in direction i Πab artificial viscosity term for interaction of particle a & b A,B,C coefficients in Johnson-Cook model ab F  Force vector on particle a due to particle b r F , z F Force component in radial and axial direction respectively r n , z n Interface normal unit vector components in radial and axial directions respectively Interface ab F  Force vector on particle a due to particle b (when a & b interact from across the interface) α, β coefficients in artificial viscosity term ab c average sound speed of particle a and b a η corrected 2D density of particle a ra radial distance of particle a from the symmetry axis a f 1 correction function for density equation of particle a a ζ non-dimensional radial distance of particle a from the symmetry axis Y σ yield stress for the material 0 Y σ initial yield stress for the material K, 1 α coefficients for the yield stress-strain curve

2 )

 2 At any time, step, as the material deforms the stress state is updated. If the von Mises stress goes beyond the Yield stress, the stress state is assumed as a trial state ( t VM σ ). The stress state then has to return back to the yield curve as shown in fig B.1. The following equation then gives the solution to the incremental plastic strain (ΔεP):

  differencing method is used to calculate the derivative of plastic strain and The above calculation can be subsisted in equation 48 to give yield stress and finally the yield stress is substituted in equation B.3 which then becomes a non-linear equation in P ε ∆ , iterative methods such as newton Raphson can then be used to solve for P ε ∆ at each time step.

Figure B. 1 .

 1 Figure B.1. Schematic shows the yield curve and the stress return algorithm in the π-plane, shows how stress from the trial state returns back to the yield curve.

Figure B. 2 .

 2 Figure B.2. Schematic shows the decomposition of elastic and plastic strain and return to the yield curve from a trial stress on a stress-strain curve.

Table 1 .

 1 Table showing material properties for A-2205.

	Parameter	Solid (A-2205)
	Density (kg/m 3 )	7800
	Sound speed (m/s)	4883.33
	Young's Modulus	186 GPa
	Poisson ratio	0.3
	A0 (in equation 15)	508 MPa
	B0 (in equation 15)	832 MPa

C0 (in equation 15) 0.031 (used as 0 for validation simulation) n (in equation 15) 0.29 0 p ε 0.05 s -1

Appendix A

Return mapping algorithm

To calculate plastic strain, a return mapping algorithm is required to return from a trial stress state to the yield curve at any time step, an incremental plastic strain is then calculated corresponding to the return from trial stress to the yield curve. Firstly, we define the von Mises stress as: