Crack nucleation in solid materials under external load-simulations with the Discrete Element Method

Numerical analysis of cracking processes require an appropriate numerical technique. Classical engineering approach to the problem has its roots in the continuum mechanics and is based mainly on the Finite Element Method. This technique allows simulations of both elastic and large deformation processes, so it is very popular in the engineering applications. However, a final effect of cracking fragmentation of an object at hand can hardly be described by this approach in a numerically efficient way since it requires a solution of a problem of nontrivial evolving in time boundary conditions. We focused our attention on the Discrete Element Method (DEM), which by definition implies ”molecular” construction of the matter. The basic idea behind DEM is to represent an investigated body as an assemblage of discrete particles interacting with each other. Breaking interaction bonds between particles induced by external forces imeditelly implies creation/evolution of boundary conditions. In this study we used the DEM approach to simulate cracking process in the three dimensional solid material under external tension. The used numerical model, although higly simplified, can be used to describe behaviour of such materials like thin films, biological tissues, metal coatings, to name a few.


Introduction
Cracking of materials is an extremely complicated process that includes processes in scales from atomic (breaking intermolecular bonds) up to a scale of thousands of kilometres in the event of catastrophic earthquakes (in the energy scale from individual eV to 10 24 J).[1] Such a large span of the scale raises a lot of questions, in particular about scalability of cracking processes, existence of factors determining the final size of the fracture area (on a macroscopic scale), course of the preceding and occurring processes during material destruction, etc [2].
The aim of this research was to try to look at the cracking processes on a scale typical for engineering and seismology (millimetres to meters) using micro-physics methods.The proposed research methodology was based on large-scale simulations using the Discrete Element technique.We mainly focused on cracking hypothetical three-dimensional materials subjected to uniaxial stretching with constant velocity and sample deformation.The above assumptions underlying the simulations may seem quite unrealistic.In fact, however, they quite well allow to describe the behaviour of cracking structures such as thin films (eg biological structures), metal coverings (eg aircraft fuselages), tailoring materials, etc.

Crack propagation
The well-known fact is that cracking solid bodies are determined by the structure of a given material, its atomic and micro-structural structure, but also by the way of applying external forces leading eventually to its destruction and fragmentation.[3] Usually, the rupture source is described by a single crack or dislocation, following the pioneering vision of Griffith [4].The dynamic crack propagation causes a relaxation of stresses and energy release, leading in the consequence to material failure.It was shown experimentally that the microdestruction leads to macrodestruction.
Correct analysis of complexity of the fracture process or/and interactions of microcracks at big concentrations typical for prefracture state is possible only in terms of statistical models.The kinetic model of evolution of crack population was introduced by Czechowski [5] and developed in [6,7,8,9,10].It lies at a level intermediate between the purely statistical approach and the fully microscopic treatment.The elementary objects are microcracks which can nucleate, propagate and coalesce.The problem of crack interaction and fusion is faced in its simplest aspects (binary interaction) but avoids its most delicate features by introducing extramechanical probabilistic assumptions.The kinetic approach operates on crack size distribution function that evolution is governed by the modified coagulation equation (mesoscopic level).Relations with the macroscopic picture, concerning the stress field evolution and the relationship between the time to fracture and the applied stress, were derived.

Discrete Element Method
Numerical analysis of cracking processes require an appropriate numerical technique.Classical engineering approach to the problem has its roots in the continuum mechanics and is based mainly on the Finite Element Method.[11] This technique allows simulations of both elastic and large deformation processes, so it is very popular in the engineering applications.In the process of cracking, the object is fragmented -disintegration into many independent bodies, and this effect is that no MATEC Web of Conferences 165, 22019 (2018) https://doi.org/10.1051/matecconf/201816522019FATIGUE 2018 numerical technique based directly on the mechanics of continuous media is clearly unable to adequately describe.Fragmentation of an object can hardly be described by this approach in a numerically efficient way since it requires a solution of a problem of nontrivial evolving in time boundary conditions.
For this reason to describe solid material fragmentations we focused our attention on the Discrete Element Method (DEM), which by definition implies "molecular" construction of the matter.[12,13] The basic idea behind DEM is to represent an investigated body as an assemblage of discrete particles interacting with each other.Breaking interaction bonds between particles induced by external forces imeditelly implies creation/evolution of boundary conditions -the most important element of the fragmentation proces.In this study we used the DEM approach to simulate cracking process in the three dimensional solid material under external tension.The used numerical model, although highly simplified can be used to describe behaviour of such materials like thin films, biological tissues, metal coatings, to name a few.
DEM is a computer simulation technique introduced by Cundall and Strack [14] and widely used in chemistry and engineering.In geotechnical application it was mostly applied to modeling granular media (soil, etc).The method is particularly convenient for rupture process study because it accommodates, in a natural way, the medium fragmentation effects.The method consists on representation of the medium by an ensemble of finite size particles (real atoms in the case of molecular dynamic simulation), which interact according to a specified interaction forces.[15,16] Fig. 1.The simplified scheme of DEM simulation.
The evolution of the system is achieved by sequential evaluating (Figure 1), for each time step, the positions, momenta and circular momenta of all particles, using Newton dynamical equation with prescribed interaction forces and assumed initial and boundary condition.[17,18] Although the approach is computationally time consuming and requires advanced HPC computational techniques, it provides the most direct insight into fully 3D dynamics of rupturing process.In this research we used the ESyS-Particle implementation of the DEM method developed at Queensland University, Australia [19] and released for the research application as the open source software.

Simulation settings
We conducted simulations assuming that the tested materials can be described as a collection of "tight" packed in a given area of spherical particles with different radii (statistical inhomogeneity of materials) interacting with each other harmonic forces that can be represented as springs connecting the means of "molecules".The simulations performed concern a relatively simple mechanism of cracking solids under the influence of uniaxial stretching.However, they showed complete compliance with experimental observations and "wider life experience".In spite of their simplicity and simplicity of the used model of interactions, they have already provided valuable information on the characteristics of the area around the tip of the developing cohesive zone.
We studied development of the kinetic model for populations of interacting different type cracks and computer simulations of models dynamics, as well as fracturing process with different regimes of generating internal stresses.
In simulations, we inserted the block of particles into the simulation object and bonded the particles together and we specified the type of interactions between bonded particles and initialise the interaction parameters.We created non-rotational elastic bonds between bonded particles as specified in the NRotBondPrms [19] parameter set: NRotBondPrms contains five parameters: a bond tag specifying which bonded particles will undergo this interaction, a unique name for the interaction group, the elastic stiffness (normalK) of the bonds, a boolean (scaling) to specify whether to scale the stiffness with particle size, and a breakDistance specifying the separation distance that must be exceeded in order to break a bond between two particles.When a particle-pair exceeds the breakDistance the bond is removed and those two particles thereafter interact according to the interactions specified for unbonded particle pairs.[19]

Maximum uniaxial stress
As it is shown in the Table 1 we prepared various samples with different elastic stiffness of the bonds between particles.Our sample (Figure 2) had 3576 particles with radii from 0.1 mm to 0.3 mm and dimensions 10mm x 10 mm x 1mm.Firstly, we conducted series of measurements with uniaxial stretching velocity equal to 1mm/s.The aim was to calculate macroscopic Young's Modulus of the samples.This is so called calibration of the DEM model.In this way we can recognize what kind of real elastic material our samples correspond to.The results are visible in the Table 1.Then, after calibration, we carried out series of measurements with three uniaxial stretching velocities: 5mm/s, 10mm/s and 50mm/s.For each one of them we calculated the value of maximum stress when main crack occurs.The results are also in Table 1.What is interesting, maximum stress weakly depends on elastic stiffness for smaller values of this parameter (500 -20000 for velocities 10 and 50, 500 -10000 for velocity 5).In this area character of cracking seems to be similar.However, for higher values of elastic stiffness there is significant increase in the maximum value of stress.At the same time, for the same values of elastic stiffness different uniaxial stretching velocities weakly influence on maximum stress.

Three types of cracking
On the basis of stress-strain curves for our results we were able to distinguish three different types of cracking: 1) type 1 when cracking is dispersed in time and space, it is impossible to extract main event or specific time of it, 2) type 2 when there is significant one big crack and some smaller 'aftercracks' and 3) type 3 with only one large crack that divides the material into two parts.These types were marked in the Table 1 and we are able to detect a certain regularity.Type 1 is always the first one, in the middle type 3 occurs and final stage is type 2. The higher the velocity the longer type 1 appears (type 3 and 2 require stronger elastic stiffness of bonds).In the Figure 3 is shown stress-strain curve typical for type 1 of cracking.It is impossible to extract single events, this is rather continuous process.Visualization is presented in the Figure 4.During uniaxial stretching a lot of cracks nucleation foci appeared in the whole area of the sample.They become bigger over time and join with each other.Second type of cracking is presented on the basis of stress-strain curve in the Figure 5.One main crack can be distinguished and some 'aftercracks' after it.In simulation (Figure 6) it is seen as one big crack which propagate quickly along sample and some smaller nucletion sources in other parts of the material which join together with main crack after some time.In the Figure 7 third type of cracking is presented on the basis of stress-strain curve.In this type of stretching one single crack appears, easily seen in the chart in the form of huge peak.In the Figure 8 one can observe that this type of crack starts from one side of the sample and propagate very fast along it.After this event the stress drops quickly to zero and sample is divided into two separated parts.

Kinetic energy
Kinetic energy provides information about total kinetic energy of all particles in every time step.This energy is connected only with linear movement of the particles because rotational ones are not considered in our simulations.On the x axis is presented strain scaled to the maximum value during simulation.Fig. 9. Kinetic energy of all particles as a function of scaled strain for uniaxial stretching with 5 mm/s velocity.
The higher elastic stiffness of the bonds the higher kinetic energy during uniaxial stretching (Figure 9).What is interesting, for the same velocity 5mm/s all curves have almost the same shape and occur for the same strain.In comparison to the Figure 9, in the Figure 10 is shown dependence for 10mm/s velocity.The energy of the cracking process is bigger, however, not two times as one would expect on the basis of velocity change from 5mm/s to 10mm/s.What is interesting and can be seen on the basis of Figures 9,10 and 11 is that strain, when peak of the kinetic energy occurs, seems not to depend on the stretching velocity and elastic stiffness of bonds.In mentioned three figures strain is scaled to maximum value, however, in real numbers, the strain is almost the same for all simulations.

Potential energy
Potential energy describes all kind of interactions connected with bonds between particles.As in case of kinetic energy this is total value in every time step.Results are shown in the Figures 12,13,14.Two interesting issues can be observed.Firstly, as in the case of kinetic energy the strain for which peak of energy occurs seems not to depend on stretching velocity and elastic stiffness of bonds.Additionally, the maximum value of potential energy strongly depends on elastic stiffness, however, does not depend very much on velocity.For instance, for elastic stiffness 200000 N/mm total maximum potential energy is equal to about 80000 J and only slightly changes when velocity increases from 5 mm/s to 50 mm/s.

Number of bonds
The initial number of bonds between particles inside the material is about 9700.The dependence between number of bonds as well as elastic stiffness and velocity is shown in the Figures 15,16,17.For velocity 5 mm/s, paradoxically, the biggest loss of bonds occurred for elastic stiffness 1000.It clearly shows that there is no clear dependence between this parameters.It is very similar for velocity 10 mm/s in the Figure 16.Quite interesting is Figure 17 where velocity 50 mm/s is shown.For this type of stretching the change in number of bonds is the most rapid, after the cracking quickly become constant.

Summary
In our research we investigated numerous parameters connected with uniaxial stretching of the sample for different elastic stiffness of the bonds between particles and velocities of stretching.We were able to distinguish three different types of cracking for materials and we calculated macroscopic Young's Modulus of the whole sample as well as maximum stress when main crack occurs.Additionally, we checked the dependence between kinetic energy of particles linear movement, potential energy of bonds between particles and number of bonds between particles as a function of stretching velocity and elastic stiffness of the bonds.
Breaking and fragmentation of solid materials is an extremely complex process, involving scales ranging from an atomic one up to thousands of kilometers in case of catastrophic earthquakes.Such a large scale span of breaking processes opens many questions concerning, for example, scaling of breaking processes, existence of factors controlling the final size of broken area, existence of precursors, dynamics of fragmentation, to name a few.
We trust that such simulations can help, soon, to understand better different processes in geophysics.There is a growing seismological evidence based on modern, high quality seismic data providing that, in spite of the large scale of failures observed in the earthquake dynamics, the breaking of rocks building the Earth's crust is controlled by processes at the micro scale.[20,21] This suggests that during rupturing there are ongoing interactions and synchronizations of processes on significantly different spatial and energy scales leading to the final post-failure state.[22] Thus, the exhaustive description of the breaking process should take these effects into account.Unfortunately, such a unified description is currently still impossible.

Fig. 2 .
Fig. 2. The sample used in the simulations which consists of 3587 randomly distributed particles with radii 0.1 mm -0.3 mm.

Fig. 10 .
Fig. 10.Kinetic energy of all particles as a function of scaled strain for uniaxial stretching with 10 mm/s velocity.

Fig. 11 .
Fig. 11.Kinetic energy of all particles as a function of scaled strain for uniaxial stretching with 50 mm/s velocity.

Fig. 12 .Fig. 13 .Fig. 14 .
Fig. 12. Potential energy of all bonds between particles as a function of scaled strain for uniaxial stretching with 5 mm/s velocity.

Fig. 15 .
Fig. 15.Number of bonds between particles as a function of scaled strain for uniaxial stretching with 5 mm/s velocity.

Fig. 16 .
Fig. 16.Number of bonds between particles as a function of scaled strain for uniaxial stretching with 10 mm/s velocity.

Fig. 17 .
Fig. 17.Number of bonds between particles as a function of scaled strain for uniaxial stretching with 50 mm/s velocity.

Table 1 .
Micro-and macro-parameters of the simulations.Maximum uniaxial stress and type of cracking are shown for different uniaxial stretching velocities and elastic stiffness of the bonds.