DEFORMATION AND FRACTURE ANALYSIS OF ELASTIC SOLIDS BASED ON A PARTICLE METHOD

In the present paper, the linear elastodynamics and fracture are simulated by using a Lagrangian particle method. The numerical model is based on the Moving Particle Semiimplicit Method (MPS) that was first developed to simulate the behavior of incompressible fluids by Koshizuka et al. (1995). The main strategy of the MPS is to replace the differential operators of the governing equations by discrete differential operators on irregular nodes, which are derived from a model of interaction between particles. In general, as a meshless method, it is very effective for the simulation of hydrodynamics problems involving free-surfaces, fragmentation and merging, and problems involving large deformation, complex shaped bodies and moving boundaries.In the last decade, the MPS method was extended to the analysis of elastic and elastic-plastic structures making possible the analysis of dynamic systems and the coupling of hydrodynamics and structure analysis to investigate hydro-elasticity problems. The implementation showed in this paper is an improved version of the simulator developed in the Numerical Offshore Tank (TPN/USP). In case of 3D analysis, instead of Euler Angles, the angles are determined by Hamilton’s quaternion algebra to avoid the singularities. On the other hand, a more generic contact searching algorithm is adopted to allow the investigation of collision and fracture amount multiple solids. The qualitative and quantitative validations of the method are carried out herein considering static and dynamic cases subjected to different boundary conditions by comparing the numerical results from MPS with Finite Element Method (FEM) and analytical solutions.


INTRODUCTION
The physics of elastic solid is based on classical elasticity theory, where the continuum field is represented by a boundary value problem and a set of partial differential equations with boundary conditions.To overcome simplified problems and solved it with exactly mathematical techniques, a set of methods have been proposed, studied and developed.The purpose of those methods is to provide an adequate discretization to domain and solve the partial differential equations.The most commons methods are the Finite Difference Method (FDM), Finite Element Method (FEM) and Boundary Element Method (BEM).Both the FEM and FDM methods are similar as the entire solution domain needs to be discretized, and a mesh is needed, while in the BEM method only the bounding surfaces need to be meshed [3].The presence of a mesh in numerical simulations can become a bottleneck in computations, therefore, meshfree and particle methods has been proposed to be computationally effective to discretize a continuum by only a set of nodal points, or particles, without mesh constraints [10].There are various meshfree and particle methods to solve elastic solids problems, such as Smoothed Particle Hydrodynamics (SPH) [11,12], Distinct Element Method (DEM) [5] and Element Free Galerkin Method (EFGM) [2].No mesh is required in meshfree particle methods, thus they are very effective for the simulation of the hydrodynamics that involve free-surfaces, fragmentation and merging, and problems that involve large deformation, complex shaped bodies and moving boundaries, in general.
The Moving Particle Semi-Implicit Method (MPS) is one of these particle methods, which was originally developed for incompressible flows with free surface by Koshizuka et al. (1995) [7].The MPS are based on Lagrangian description and the computational domain is discretized by particles that interact among themselves and they according to the governing equations of interest.For elastic solids, the constitutive equation characterized by the material density, Young's modulus and Poisson's ratio are replaced by algebraic equations derived from a model of interaction between particles.This results, essentially, a system of particles connected by normal and tangential springs.The motions of translation and rotation of each particle are calculated by numerical integration of the equation of motion.
The implementation showed in this paper is an improved version of the simulator developed in the Numerical Offshore Tank (TPN/USP), based on studies [4,8,13,14].In case of 3D analysis, instead of Euler Angles, the angles are determined by Hamilton's quaternion algebra to avoid the singularities.On the other hand, a more generic contact searching algorithm is proposed to allow the investigation of collision and fracture amount multiple solids.To validate the method, the qualitative and quantitative simulations are carried out herein considering static and dynamic cases subjected to different boundary conditions and the comparison of numerical results from MPS with FEM and analytical solutions.

Governing equations
The governing equations of the dynamic of elastic solid can be written as [14] where r α is the position vector, v α is the velocity vector, ρ is density, ε αβ and ε γγ are components of strain tensor, δ αβ is the Kronecker's delta, b α is the body force vector and α, β and γ are the components dimensions.The Lamé's constants µ and λ are given by Eqs. ( 3) and (4), respectively.
where E is the Young's modulus and ν is the Poisson's ratio.The equation of motion Eq. ( 2), can be rewritten by introducing stress tensor σ αβ and isotropic pressure p as Here,

Numerical model
In MPS method, the differential operators of the governing equations are replaced by discrete differential operators on irregular nodes [6].For a given particle i, the influence of a neighbor particle j is defined by weight function w(| r ij |) given in Eq. (8).
where r e is the effective radius that limits the range of influence and r ij is the distance between i and j.For the cases simulated in this paper the value of r e is 2.1l 0 , where l 0 is the interval between the adjacent two particles in the initial configuration.As a result, for a function φ, the gradient, divergent and rotational operators are defined in Eqs. ( 9), ( 10) and (11), respectively.
where d is the number of spatial dimensions, s ij is a versor perpendicular to r ij and n i is the particle number density defined as

Interaction between particles
Initially two particles are separated by r 0 ij .After the motion, the new position vector is represented as where x i and x j are the positions of the particles i and j.
The displacement vector between i e j can be written as where R ij is the rotation matrix, in 2D written as and in 3D using quaternion q as [1], we have where (η x η y η z ) is the versor η ij representing the axis of a rotation and defined by the average angular velocity vector ω ij , Eqs. ( 19) and (20).
The quaternion at the next time step q k+1 is obtained as The displacement vector can be divided into normal u n ij and shear u s ij components, Eqs. ( 22) and ( 23) , and the strain vector is calculated as Eqs.( 24) and (25).
The volumetric deformation ε γγ is described using the divergent of normal displacement vector.
Figure 1 shows the displacements and rotations between particles.
The normal stress vector σ n ij , shear stress vector σ s ij and isotropic pressure p i are obtained according to the deformation components.
Using the divergent operator to Eq. ( 5), translation of particles is obtained from normal stress vector, shear stress vector and isotropic pressure.
where p ij is the average pressure, Eq. (33).
The shear stress vector also influence in the rotation of particles, where the force is calculated as Eq. ( 34) and the moment can be written as Eq. ( 35).
If the moment of inertia I i is constant along the time, the angular acceleration vector of particles is calculated as

Velocity and position
New velocity v k+1 i and position r k+1 i are calculated by using a explicit Euler method, Eqs.(37) and (38).
Also by the explicit Euler method, new angular velocity w k+1 i and rotation θ k+1 i can be calculated as where ∆t is the time step.

Fracture
In order to simulate multi-body dynamics with fracture, a more generic algorithm for the detection of contact is proposed herein.
Initially the particles of solids are divided into surface particles and internal particles.Similar the condition of free surface for fluid [9], the particle number density n i is used as criterion.If the value of n i is lower than a certain value, Eq. (41), the particle is set to surface particle.
where n max i is the maximum value of particle number density and n c is a constant value.In the present work is used n c = 0.5.
When the solids move closer each other, the surface particles detected different solids and the range of influence r e is set equal to initial particle distance l 0 .With the range of influence delimited by initial particle distance, the interaction between surface particles occurs only when they are in compression and thus only repulsive forces acting between surface particles to avoid the overlapping of particles.In case of fracture, the occurrence of fracture is detected when the strain ε ij between particles i and j is greater than a critical distance of fracture ε max .In this situation, the weight function is set to zero for the particles i and j subjected to fracture, Eq. ( 42), and they are set to surface particles.This will make the particles i and j ready for detection of contact for the next time step.Figure 2 gives an example of fracture of an elastic body and the change of internal particles to surface particles during fracture.

RESULTS/VALIDATION
Three cases are simulated and compared to verify the model implemented.First a static case, obtained by gravity imposed gradually, is simulated and the displacements of a fixed-free beam are compared with analytical results.A dynamic case involving the motion of a fixed-supported beam is simulated and the maximum deflection, the natural frequency and the period of vibration obtained from MPS are compared to analytical and numerical results.Finally a collision case shows fracture with contact detection proposed in this paper.

Static case
Figure 3 shows a fixed-free beam of length 1.0 m and square cross section of 0.2 x 0.2 m subjected to self-weight along longitudinal axis.The transversal displacements dx and dy are given by Eqs. ( 43) and (44), respectively, and the longitudinal displacement is given by Eq. (45).The material properties are density ρ = 1000 kg/m 3 , Young's modulus To compare the results from MPS, which is implemented for dynamic analysis, to analytical results, wich is applied to static cases, a function given by Eq. ( 46) is imposed gradually simulating the gravity acting in the beam, as shown in figure (4).Thus, the comparison can be made when the steady state is reached.
where G is the gravity magnitude = 9.81 m/s 2 .Figure 5 shows the results of the displacements dx and dz in x and z directions, respectively.The displacement dx is measured at x = 0.1 m of the cross section and dz is measured at x = 0.1, y = 0.0 m of the cross section.Horizontal axis of figure 5 is the position z of the cross section and vertical axis is the displacement.The results shows that the computed dx agrees very well with the analytical ones.The error for dz is slightly larger, but in general agrees well with analytical solution.

Dynamic case
As a simple dynamic case, the first mode shape of a fixed-supported beam of length 1.0 m with square cross section of 0.225 x 0.225 m, as shown in figure 6, is considered.The results from MPS are compared to FEM and analytical results.The first mode shape is excited with gravity and the motion amplitudes is given in Eq. (47).
The plotted results of deflection can be seen in the graphic presented in figure 7. Table 1 shows the results of natural period, frequency and amplitude for the first mode shape.The natural period obtained from MPS presents a small discrepancy if compared to the result from FEM and a considerable difference if compared to analytical result, where the relative error is about 11.11%.The maximum deflection is closer to the analytical result with relative error about 10.71%, while the relative error between MPS and analytical is slightly larger, about 19.02%.The results show good agreement despite deviations that can be considered as numeric errors and simplifying assumptions.The collision of two elastic solids occur at 0.07 s.After the fracture occurs in the left side of the block at 0.60 s, the collision between particles of the fractured surfaces is detected at 0.70 s.The repulsion due to collision between particles of fractured surface is visible in the case with contact detection at 1.00 s.As a result, the parts of block move away due the collision at 0.70 s.On the other hand, as shown in the left column of figure 9, if detection between the surface created by fracture is not detected, a new and stronger collision occurs at 1.70 s and the simulation stopped owing to the overlap of the particles.

CONCLUSION
A computer code to simulate the dynamics of elastic solids has been implemented.Simulations are carried out for classic examples of beams and the results from MPS are consistent with the numerical and analytical results.A detection condition is proposed to avoid overlapping of particles due the contact between the fractured surfaces.A simulation of collision involving fracture shows problems if the detection of particles is neglected in fracture case.

E = 10 M
P a and Poisson's ratio ν = 0.3.The simulation parameters are particle distance dp = 0.025 m and time step ∆t = 10 −5 s.

Figure 8
Figure8shows initial conditions of the case of collision between two elastic solids involving fracture.Figure9gives the snapshots of the simulation.The left column shows fracture without contact detection and the right column shows fracture with contact detection.A cube (salmon) 0.5 x 0.5 x 0.5 m with an initial velocity v y = −5 m/s collide into a block (blue) 1.0 x 0.25 x 2.0 m initially stopped.The material properties of cube are density ρ = 1000 kg/m 3 , Young's modulus E = 6 M P a and Poisson's ratio ν = 0.3 and the material properties of block are density ρ = 1000 kg/m 3 , Young's modulus E = 10 M P a and Poisson's ratio ν = 0.3.The simulation parameters are particle distance dp = 0.05 m and time step ∆t = 10 −6 s.The critical distance of fracture is ε max = 0.2.

Figure 8 .
Figure 8.Initial conditions of the case of collision.

Figure 9 .
Figure 9. Left: without contact detection, Right: with contact detection.

Table 1 .
First mode results