A unified approach for static and dynamic fracture failure in solids and granular materials by a particle method

The material structure at the microscale reveals the particulate nature of solids. By exploiting the discrete aspect of materials, the so-called particle methods developed and applied to the simulation of solids and liquids have attracted the attention of several researchers in the field of computational mechanics. In the present paper, a particle method based on a suitable force potential is proposed to describe the nature and intensity of the forces existing between particles of either the same solid or different colliding solids. The formulation applies to problems involving both granular materials and solids interacting with granular materials. The above approach is applied to simulate different problems dealing with the 3D dynamic fracture and failure of solids.


INTRODUCTION
he particle methods have attracted the attention of several researchers in the field of computational mechanics [1,2].It is straightforward when the material has to be described at the nanoscale where the molecular dynamic method can easily be applied [3]; however, this discrete description can successfully be applied even at the macroscale [4,5].Accordingly, a discrete approach can replace the traditional continuous modelling of a solid, and allows us to examine complex phenomena such as failure, fracture, and interaction with other bodies [4].The formulation applies to problems involving both granular materials [6][7][8] and solids interacting with granular materials.The need to solve mechanical problems has solicited the development of several numerical techniques aimed at assessing the response of materials both during service and at their ultimate conditions.To this purpose, various computational tools have been proposed: the finite element method (FEM), the finite difference method, the finite volume method, the more recent boundary element method, meshless method [9] and natural element method [10], just to cite the main tools.Problems involving mechanical systems characterized by large displacements, high velocity impact, fracture, failure, fragmentation, clustering, fluid-solid interaction, and so on, cannot be solved through classical numerical techniques in the context of the so-called Lagrangian approach, because the severe distortion of the discretized domain leads to a lack of consistency between the numerical modelling and the physics problem.This issue has partially been overcome by using special remeshing techniques and the FE method, but that presents some drawbacks.The discretization through only nodal points -i.e.without mesh connectivity -represents a possible approach to solve this problem.In this context, the so-called Smoothed Particle Hydrodynamics (SPH) has been one of the first particle meshless methods in computational mechanics [11,12].The original idea of the developers of the above method lies in T the fact that the collective motion of particles in particle-wise systems (such as the movement of liquid, gas flow, particulate assemblies, etc.) may physically be modeled by the governing equations of classical Newtonian hydrodynamics.Such mesh-free Lagrangian methods are characterized by variable nodal connectivity; the material is discretized with interacting particles and the method transforms a discrete field into a continuum one by using a localized kernel function that operates as a smoothing interpolation.The method can be applied to discrete systems but also to continuum ones, where the nodal points are not directly associated with the single particle composing the matter [13][14][15].The noticeable advantages of such a method became clear after its first introduction, suddenly leading to its extension to the solution of a wide range of physics and engineering problems, well beyond the initial classical hydrodynamics field [16,16].The discrete nature of solid matters is evident at the microscale or at the molecular level [17].On the other hand, the matter is typically modeled as continuous medium at the macroscale, suggesting the use of continuous computational techniques such as the FEM.However, discrete models can conveniently be applied even in macroscopically continuum systems where each representative particle can be associated to a cluster of continuum elements.Discrete Element Methods (DEMs) sometimes also called Particle Methods (PMs), are numerical procedures to solve different engineering problems: the main assumption for this class of methods consists in the description of the material through an assemblage of separate discrete elements, such as atoms, molecules, grains, particle solid elements, etc. according to the scale of observation adopted [1,17].In the discrete approach, non-linear interactions between bodies and within bodies are numerically simulated in order to get their non-linear differential equations of motion.It is worth mentioning that the discrete nature of materials is well evident for some classes of them, such as the granular ones: geomechanical and powders materials can be studied by exploiting their noticeable granular nature [18,19].Discrete methods have usefully been applied to solve problems at the macroscale such as mineral processing, rock and concrete blasting and crushing, sand mechanics, failure and fracture of compact or granular bodies [20,21], problems involving fluid materials like liquids and gases [16,22].A generic solid can also be assumed to have a particle structure with a proper nature of particles' interaction forces.From this remark, the mechanical behaviour of a material can range from very incoherent behaviour (characteristic of fluids, granular, powders, etc) up to compact behaviour (typical of polycrystalline materials), respectively.Therefore, a discrete model can be used to simulate different classes of solids, by properly choosing the nature of the interaction forces between their discrete constituent elements.Such an approach permits to get the overall response of the material at the macroscale, which is the main interest of materials science and mechanics of materials.Moreover, the particle discretisation for a solid allows us to tackle the problem from a dynamic point of view, enabling to solve high strain rate, impact, large displacements and large strains problems involving history-dependent behaviour, plasticity, etc., that are more naturally expressed in a Lagrangian computational framework.It should also be considered that discrete methods allow us to deal with short range forces (such as contact forces) and long range forces (such as the electrostatic ones).That enables to represent the physical nature of materials at the small or nano scales, where interatomic forces obeying such a non-local interaction relationship exist.In the present paper, a discrete element approach for the mechanical simulation of either continuum or granular-like materials (formulated on the basis of the force potential interaction concept) is proposed and applied to analyse dynamic fracture and failure of both granular materials and compact solids.A brief introduction related to the forces existing in continuum solids and discrete incoherent aggregates -modeled through discrete particles -is presented by taking into account the dynamic nature and large strain characteristics of the problem.Finally, some examples demonstrating the versatility and the capability of the proposed approach are discussed.

POTENTIAL-FORCE FORMULATION FOR PARTICLE-PARTICLE INTERACTION
s is recalled above, the particle approach is straightforward when the material is really composed by several discrete elements as at the nano/microscale (for both solids and granular-like matters), where the particles can immediately be identified as the atoms/molecules or grains in the matter.However, a continuous solid can also be idealized as a cluster of elements reciprocally joined by forces depending on the particles relative positions: if the particles tend to approach to each other, a repulsive force arises (typical of compact solids and granular materials), while an attractive force appears when they tend to move far away (typical of compact solids and cohesive granular materials).The above interacting forces can be thought to have a local nature when they are associated to the direct contact between elements (for instance when granular materials are examined), while they assume a longer range characteristic in the case A of continuous materials.This remark suggests to adopt a unified potential-based formulation for the assessment of the interparticle actions, by properly choosing the distance of influence for the transmission of forces.As is typically done in the atomic description of solids where the material can be characterized by a potential energy functional ( )  x , it can be assumed the existence of a similar potential function at the meso or macro-scale.The equilibrium equations, referred to the discrete element i, correspond to the configuration of minimum energy of the system: ( ) 0 where ( ) tot  x and i P are the strain energy and the force applied to the particle i, respectively.By writing the power series expansion of the potential energy function, starting from its equilibrium state identified by the position vector 0 x , the stationary condition of the above functional identifies the equilibrium configuration as follows: where x x 0 at the equilibrium, and the stiffness matrix of the system is indicated with At the atomic scale, various potential functions have been proposed in order to represent the mechanical behaviour of materials, such as the Lennard-Jones (LJ) potential and the Morse potential [23].The simplest mechanical model describing the stress-strain behaviour of a continuum linear elastic solid is the generalized Hooke's law that assumes the force between two infinitely closed material points to be proportional, through proper coefficients, to the strain value estimated in such a small region (Fig. 1).The region of material -representative of a small region (of a continuum solid) lying between two particles i, j -is subjected to the deformation (evaluated with respect to the equilibrium distance e r ) given by:  Eq. ( 3) is written by taking into account the large displacements effect, whereas the maximum co-penetration distance enables to avoid the complete penetration condition between particles.By assuming the linear stress-strain relationship, the following expression for the force acting along the line joining the particle centers can be written: where ij A and 0 '( ) ( ) E s c s E   are the cross-section and the elastic modulus of a bar element that represents the material connecting the particles , i j .The Young modulus '( ) E s is assumed to depend on s in order to take into account the no-penetration condition: it increases when 0 s  (i.e. ( ) c s   if 0 s  ), whereas ( ) 1 c s  when s increases.The assumption of the no-penetration condition ( 0 ) entails an unlimited compressive strength of the particles when they are in contact, up to the limit 0 s  (i.e. 2 w   or 0 r r  ), w being the particles overlapping and 0 r being the minimum allowed distance between particles.The force-potential and the corresponding stiffness, according to Eq. ( 3), can be written as follows: e ij e e e s s s s A E s s s s In Eq. ( 3) the influence radius infl r indicates the maximum interacting distance between two particles, i.e. it is the maximum distance beyond which two particles are not interacting.This assumption is suitable to represent long distance forces, such as in continuum materials, since the internal actions between the material points act also if such points are not in direct reciprocal contact (non-local interaction).In Fig. 3a the distance dependence of the potential, of the force and of the stiffness are displayed for the linear elastic particles interaction related to the representative equilibrium distance / 1 e i s d  , while the same functions are reported in Fig. 3b for the case of no-penetration constraint.In order to correctly represent the elastic behaviour of the material described through bar elements connecting the particles, a suitable choice of the cross section area ij A of such elements must be done.Since the bar arrangement in the 3D space depends on the influence radius infl r and on the particles distribution (the bar connecting two particles exists only if the particles are within their influence distance), we can assume infl 0, ( ) , where the reference area is The above potential formulation can also be adopted for granular materials, where the internal forces act only when the particles are in reciprocal contact and their motion make them approach.In these cases, an influence distance can be adopted ( infl / ( / 2) i r d lower than or equal to about 2), i.e. the interaction of particles occurs only when they are in contact under compression condition, while no forces exist when ' 0 s  .Further, a tangential force can be assumed to exist between the colliding particles when ' 0 s  (Ref. [22]):   , , ( ) sign( ) min , ( ) where  is the viscosity coefficient, md  is the dynamic friction coefficient of the materials, , t rel u is the relative tangential velocity between the two particles in contact.Note that Eq. ( 5) provides a tangential force as the minimum value between the dynamic friction and the viscosity action.The present authors have verified that the function infl ( ) a r , for a given particles arrangement, is almost independent of the particles arrangement in the space and depends on infl r only [24].On the other hand, the Poisson's coefficient tends to 1/3 as the influence radius increases, irrespective of the particles layout used to represent the domain of the problem of interest.When compact materials have to be simulated, a suitable choice of the influence radius value is infl / ( / 2) 4 for both cubic and tetrahedral arrangements), i.e. the interaction forces exist for neighbor particles even if they are not in direct contact.Granular materials can be simulated by assuming infl / ( / 2) 2 i r d  and infl ( ) 1.0 a r  , i.e. the interaction of particles exists only when they are in reciprocal compressive contact, while no forces are present when ' 0 s  .It is worth noting that the particle-boundary simulation can be tackled in the same way as the particle-particle interaction.

SPH property of the potential force approach
Smoothing Particle Hydrodynamic (SPH) approximation of a function and of its n-th derivative at a point i x in the region  can be written as follows [13]: where ( ) i   x x is the Dirac delta function that is numerically approximated by adopting a proper decreasing function ( ) i W  x x of the distance i  x x (W is also called kernel or smoothing function), and h is the smoothing length (or support dimension).The numerical evaluation of Eq. ( 61) is expressed as follows: where the index k spans all the particles over the neighbors of that placed at i x because of the compact support of the kernel function.The parameters , , x are the position, mass and density of the generic particle k , respectively, and ( ) k f x is the function value at k x .By using a SPH approximation for the material's stress field and inserting it in the dynamic equilibrium equation, , 0 , the following expression can be obtained: where p V represents the volume of the particle p and / . By using the potential formulation presented above, we get: where ( ) i kp q is the i-th component of the versor connecting the particles k and p, and [ ] k p indicates the generic particle k falling within the neighbor of particle p.

DYNAMIC EQUILIBRIUM EQUATIONS
he dynamic governing equations of the discretized problem are given by: where M is the mass matrix, x  is the vector of the particles center acceleration (the rotation degrees of freedom are neglected since particles with small rotational inertia are assumed).Further, , , F are the internal force vector (see previous sections), the damping force vector, the external force vector and the vector of actions due to the collision of the particle with the elastic boundaries, respectively.The numerical integration of the motion equations can explicitly be done for each particle, once the force vector acting on it is known.Explicit time integration methods, even if conditionally stable, are usually adopted in particles simulations to avoid the use of large matrices such as the stiffness one.T In the present particle approach, the explicit leapfrog method (the Verlet integration scheme) is adopted [25], which allows a satisfactory numerical stability, provides time-reversibility and preservation of the symplectic form on the phase space.The word leapfrog recalls that even derivatives of position are known at on-step points, whereas odd derivatives are known at mid-step points.The particle acceleration at time step n , the velocity at time step 1/ 2 n  and the position of particle i at time step 1 n  are given by: , ,, / is the mean velocity vector across the time step 1 n n   .In order to determine stable results, an upper limit of the time step size / t m K     in the DEM approach is generally obtained from the natural oscillating period of a couple of particles: in the above expression, the minimum ratio / m K among all the particles constituting the discretized solid must be taken into account, and 1   is an appropriate constant to be introduced.

Granular material cutting simulations
his example considers the cutting of a granular material that resembles the fracture phenomenon in solids.A vertical flat blade with a constant horizontal velocity equal to 10 mm s -1 [26] drives particles that move, leading to some failure lines in the body.An initial volume of particles (Fig. 4) under the gravity action is assumed to be placed inside a box, and the vertical blade is placed at the left side position.The simulation consists in moving the blade, causing the material inside the box to redistribute.The performed dynamic analysis is aimed at determining the particles configurations at different time instants.The material of the particles (assumed with a mass density equal to 855 kgm -3 ) is supposed to have an elastic modulus equal to , and the two coefficients of dynamic friction (between particle and particle, and between particle and planes) are supposed to be 0.10 . The obtained arrangement of the particles shows a behaviour similar to that of the soil cutting case analysed in Ref. [26].Note that the blade in the present study is shorter than that in the literature test, causing a little difference in the disposition of the particles near the blade.Neglecting last aspect, the shape assumed by the particles representing the granular material in the present study is similar to the configuration presented in Ref. [26].In Fig. 4e, f, the fracture lines arising in the material can be recognized, i.e. the cutting phenomenon can be thought as a crack appearance event inside a non-compact material.

Fracture of a brittle concrete beam under a dynamic impact
This example examines the fracture behaviour of a compact material.This analysis is performed with the same potential approach described in the theoretical Sections and used also for the previous numerical example.In the present numerical test, a simply supported plain concrete beam under impact load is taken into account [27]  , mass density equal to 8000 kgm -3 ).The system is modelled with about 2200 particles with cubic arrangement and radius equal to 35mm each.The impact velocity of the steel mass is assumed to be equal to  Fig. 5 shows the configuration of the system at different time instants.In particular, the initial development of the elastic wave -propagating inside the beam -is shown in Fig. 5a, whereas the failure pattern at t=0.125 s after the first impact is represented in Fig. 5b.As can be observed, a wedge-like failure mechanism takes place producing the detachment of a wide bottom portion of the beam.Moreover, a diffused separation of the bottom layer of the beam can also be acknowledged.In Fig. 5c, the time history of the reaction force at support obtained by the discrete approach by using different discretizations is displayed and compared with the FE results provided in Ref. [27].It can be observed that, only after a first time interval equal to about 30 ms (when the reaction is practically zero due to the time required by the elastic wave to propagate from the stricken zone to the support), such a boundary force raises very quickly.Of course, higher impact speeds correspond to higher reaction force values.The fracture pattern provided by the DEM approach roughly corresponds to the literature results, despite the other FEM approach is very different (Fig. 5d) and does not allow us to determine the detachment of material portion produced by the fracture phenomenon.

CONCLUSIONS
he particulate nature of solids, naturally acknowledged at the microscale in granular or fluid materials, can also be adopted conveniently for the mechanical description of compact matter at the macroscale.Based on such an idea, a particle method can be applied to the simulation of a wide class of mechanical problems once the mechanical properties and governing equations have been established.A discrete approach can replace the traditional continuous modeling of solids, allowing to examine complex phenomena such as failure, fracture, contact and interaction with other bodies, by also taking into account large strain and dynamic effects.In the present study, a particle method -based on the description of the forces between discrete elements evaluated on the basis of a force potential -has been proposed.It has been emphasized that the formulation applies without changes to problems involving granular materials or solids interacting with granular materials.The above approach has been applied to the simulation of different problems dealing with the 3D fracture and failure of both granular materials and compact solids.The formulated discrete method shows a wide capability to deal with different complex mechanical problems by properly describing the evolution of the crack pattern under dynamic conditions.
the distances between the particle centers in a generic state and at the equilibrium state, respectively, whereas , e s s are the corresponding effective distances counterparts (Fig.2.In order to account for the no-penetration condition, the maximum co-penetration depth 1 2 r         is expressed as a fraction  of the averaged radius r of the particles.

Figure 2 :
Figure 2: (a) A couple of particle at the equilibrium distance.(b) Maximum co-penetration corresponding to

Figure 3 :
Figure 3: Force potential and corresponding forces and stiffness according to the linear elastic behaviour for an equilibrium distance / 1 e i s d  , (a) without and (b) with the no-penetration condition.

8 3
tensile strength (or, equivalently, no cohesion).The force potential is used to describe the particles interaction by adopting an influence distance of the particles equal to infl r d  , and the box and the blade are supposed to have elastic modulus equal to volume is modeled through about 1162 spheres with diameter size equal to 20 mm, initially arranged in a cubic fashion.The geometry of the system is depicted in Fig.4.The integration time step is taken equal to 0.2 ms t   and the final blade displacement is equal to 200mm, with an analysis duration equal to 2.0 s.The damping coefficient for the dynamic analyses is assumed equal to 0.20d   equal to 2300 kgm -3 for the concrete material, whereas the falling hammer is supposed to be made of steel (with elastic modulus potential described above is used to quantify the particles interaction, by adopting an influence distance of the particles equal to infl 2 r d  and a time increment for the dynamic analysis equal to 5

Figure 4 :
Figure 4: Configuration of particles for a blade displacement equal to: (a, c) 10cm, (b, d) 20cm [25], (e, f) corresponding configurations provided by the present model.Color scale indicates the horizontal displacement of the particles expressed in m.

Figure 5 :
Figure 5: (a) Elastic wave propagation and (b) failure pattern for the case 0 2 / v m s  .(c) Time history of the reaction force at the beginning of the impact phenomenon for 0 2, 4, 6, 8 / v m s  from the present model.(d) FE crack path reported in [26].