Numerical modeling of rock fracture and fragmentation under impact loading using discrete element method

Abstract

The fracture and fragmentation of rock materials are basic and important problem in geomechanics and blasting engineering. An approach, which can simulate the process of fracture and fragmentation of rock materials, is introduced in this work. A beam–particle model is first introduced in the frame of the discrete element method. In the beam–particle model, the neighboring elements are connected by beams. Consequently, a beam network is formed in the particle system. The strength characteristics of rock materials are reflected by the beam network. The strength criterion was then built to verify whether a beam exists or not. The process of rock fracture and fragmentation is described by the gradual disappearance of beams. Finally, two cases were presented to indicate the validity of the method proposed in this work.

Introduction

The fracture and fragmentation of rock materials are a current subject of keen interest, which are widely used in sciences and engineering, for example, mining exploitation, tunnel piercing, and protection works. At present, the numerical methods, which were used to simulate the fracture and fragmentation of rock materials, can be divided into two types: the method based on the continuum mechanics and the method based on multibody dynamics.

For the method based on the continuum mechanics, most simulations were completed by the finite element method (FEM) in which a constitutive relation and strength criterion need to be selected or constructed.^1–3 However, the construction of the constitutive relation and strength criterion is difficult since the rock materials exhibit typically discontinuous, anisotropic, and heterogeneous characteristics. Therefore, the difference exists between the simulated and the experimental results. Otherwise, the flyrock, crack formation, and development are also difficult to be simulated by the method based on the continuum mechanics. The reason lies that the common nodes between the neighboring elements in the unmixed FEM would not be separated. Nevertheless, some researchers still tried to develop different types of mechanical models in the frame of FEM to break through the limitations of FEM such as the fictitious crack model (FCM), ⁴ the discrete cracking model (DCM), ⁵ the smeared cracking model (SCM),^6,7 and the random lattice model (RLM). ⁸ The above-mentioned models are nonlinear fracture mechanics representation, which are general model to all kinds of materials. However, practical and numerical difficulties prevent its application in practical cases.

For the method based on multibody dynamics, the discrete element method (DEM) ⁹ was employed to simulate the fracture and fragmentation of rock materials since the fracture and fragmentation are the process from integrated state to granular state. The DEM, also called a distinct element method, is any of a family of numerical methods for computing the motion and effect of a large number of small particles. Although DEM is very closely related to molecular dynamics, the method is generally distinguished by its inclusion of rotational degrees of freedom as well as stateful contact and often complicated geometries. With advances in computing power and numerical algorithms for nearest neighbor sorting, it has become possible to numerically simulate millions of particles on a single processor. At present, DEM is becoming widely accepted as an effective method of addressing engineering problems in granular and discontinuous materials, especially in granular flows, powder mechanics, and rock mechanics. The DEM can realistically simulate the contacts, crush, and collisions between bodies.

A beam–particle model (BPM), which was proposed by Kun and Herrmann, ¹⁰ was employed to simulate the fracture and fragmentation of rock materials. In this work, the BPM proposed by Kun and Herrmann was extended to simulate three-dimensional problem. A beam network was first constructed to reflect the strength of the characteristics of rock materials. Then, the mechanical behaviors of beam network were described by the matrix displacement method (MDM).

Contact forces and particle motions

The simplified expressions of Mindlin’s theory^11,12 are introduced to calculate the normal and tangential forces for two particles in contact with this study. The motions of particles are described by DEM. The equations of interactions and motions are listed in Table 1.^13,14

Table 1.

Equations of interactions and motions for particles in DEM.

Variables	Equations
Normal forces	${\mathit{F}}_n^{\mathrm{ij}}=-[\frac{4}{3}·\frac{E_i{E}_j}{E_i(1-v_i^2)+E_j(1-v_j^2)}{R}_{\mathrm{ij}}^{1/2}{\delta }_{\mathrm{ij}}^{3/2}+{\eta }_d{\mathit{V}}_{\mathrm{ij}}·{\mathit{n}}_{\mathrm{ij}}]{\mathit{n}}_{\mathrm{ij}}$
Tangential forces	${\mathit{F}}_t^{\mathrm{ij}}=-[\frac{16}{3}·\frac{G_i{G}_j\sqrt{R_{\mathrm{ij}}{\delta }_{\mathrm{ij}}}}{G_i(2-v_j)+G_j(2-v_i)}{\tau }_{\mathrm{ij}}+{\eta }_{\tau }{\mathit{V}}_{\mathrm{ij}}·{\mathit{t}}_{\mathrm{ij}}]{\mathit{t}}_{\mathrm{ij}},\left|{\mathit{F}}_t^{\mathrm{ij}}\right|=min(\left|{\mathit{F}}_t^{\mathrm{ij}}\right|,\mu \left|{\mathit{F}}_n^{\mathrm{ij}}\right|)$
Translational motion	$m_i{\stackrel{··}{\mathit{X}}}_i=\sum _{j=1}^N({\mathit{F}}_n^{\mathrm{ij}}+{\mathit{F}}_t^{\mathrm{ij}})+m_i\mathit{g}$
Rotation motion	$I_i{\stackrel{··}{\mathbf{\theta }}}_i=\sum _{j=1}^N({\mathit{d}}_{\mathrm{ij}}\times {\mathit{F}}_t^{\mathrm{ij}})$
Relative velocity	${\mathit{V}}_{\mathrm{ij}}={\stackrel{·}{\mathit{X}}}_i-{\stackrel{·}{\mathit{X}}}_j+{\stackrel{·}{\mathbf{\theta }}}_i\times \left(r_i{\mathit{n}}_{\mathrm{ij}}\right)-{\stackrel{·}{\mathbf{\theta }}}_j\times \left(r_j{\mathit{n}}_{\mathrm{ij}}\right)$
Unit normal vector	${\mathit{n}}_{\mathrm{ij}}=\frac{X_i-X_j}{|X_i-X_j|}$
Unit tangential vector	${\mathit{t}}_{\mathrm{ij}}=\frac{{\mathit{V}}_{\mathrm{ij}}-({\mathit{V}}_{\mathrm{ij}}·{\mathit{n}}_{\mathrm{ij}}){\mathit{n}}_{\mathrm{ij}}}{|{\mathit{V}}_{\mathrm{ij}}-({\mathit{V}}_{\mathrm{ij}}·{\mathit{n}}_{\mathrm{ij}}){\mathit{n}}_{\mathrm{ij}}|}$
Normal overlap	${\delta }_{\mathrm{ij}}=r_i+r_j-({\mathit{X}}_i-{\mathit{X}}_j)·{\mathit{n}}_{\mathrm{ij}}$
Tangential overlap	${\tau }_{\mathrm{ij}}=({\mathit{V}}_{\mathrm{ij}}·{\mathit{t}}_{\mathrm{ij}})\Delta t$

$r_i$ and $r_j$ , $E_i$ and $E_j$ , $v_i$ and $v_j$ , $G_i$ and $G_j$ , ${\stackrel{·}{\mathit{X}}}_i$ and ${\mathit{X}}_j$ , ${\stackrel{·}{\mathit{X}}}_i$ and ${\stackrel{·}{\mathit{X}}}_j$ , and ${\stackrel{·}{\mathbf{\theta }}}_i$ and ${\stackrel{·}{\mathbf{\theta }}}_j$ are, respectively, the radius, elastic moduli, Poisson’s ratios, shear moduli, position vectors, velocities, and angular velocities of particles $i$ and $j$ ; ${\eta }_d$ and ${\eta }_{\tau }$ are the normal and tangential damping coefficients, respectively; $m_i$ and $I_i$ are correspondingly the mass and moment of inertia of particle $i$ ; ${\mathit{d}}_{\mathrm{ij}}$ is a vector from the mass center of particle $i$ to the contact point; $\mathit{g}$ is gravitational acceleration; $N$ is the number of particles in contact with particle $i$ ; $\mu$ is the friction coefficient; ${\stackrel{··}{\mathit{X}}}_i$ and ${\stackrel{··}{\mathbf{\theta }}}_i$ are the acceleration and angular acceleration of particle $i$ ; $R_{\mathrm{ij}}=r_i{r}_j/(r_i+r_j)$ ; and $\Delta t$ is the time step.

BPM

The simulated object is discretized into an assembly of particles in the BPM. In this work, the particle with spherical shape is used to discretize an object. The contact detection method ¹⁵ is first initiated to identify all neighboring particles for a particle. Then, a beam is added to connect arbitrarily two neighboring particles. As a result, a beam network is formed in the particle system. The strength of the material is reflected by the beam networks in BPM. The length of a beam is equal to the distance between two particles’ centers which lie at the beam ends. The beam cross section is supposed to be circular. The diameter of a beam cross section is determined as one-sixth of the beam length to ensure the requirement by Euler–Bernoulli beam theory. The interactions and motions of particles can be calculated by DEM. ⁹ Figure 1 shows a beam network in a particle system, and Figure 2 shows a deformed beam under the forces applied on the two ends of the beam.

Figure 1.

Beam network in a particle system.

Figure 2.

Schematic diagram of a deformed beam.

In the BPM, the deformation of a beam between two neighboring particles at a certain time is already known according to the DEM. Therefore, the end forces of a beam can be calculated by the MDM. ¹⁶ The forces and torque at the two ends of a beam can be expressed as

$$\begin{array}{cc}\multicolumn{1}{c}{}& \left\{F_m\right\}=\hfill \\ \multicolumn{1}{c}{}& {\left[\begin{array}{cccccccccccc}\hfill X_i^m\hfill & \hfill Y_i^m\hfill & \hfill Z_i^m\hfill & \hfill M_{x,i}^m\hfill & \hfill M_{y,i}^m\hfill & \hfill M_{z,i}^m\hfill & \hfill X_j^m\hfill & \hfill Y_j^m\hfill & \hfill Z_j^m\hfill & \hfill M_{x,j}^m\hfill & \hfill M_{y,j}^m\hfill & \hfill M_{z,j}^m\hfill \end{array}\right]}^T\hfill \end{array}$$

(1)

where $i$ and $j$ are two ends of the beam, respectively; $m$ is the label of the beam; $X$ , $Y$ , and $Z$ are the forces along the three coordinate axes, respectively (as shown in Figure 2); and $M$ is the torque on both ends of the beam. The relative deformation can be expressed as

$$\begin{array}{cc}\multicolumn{1}{c}{}& \left\{{\delta }_m\right\}=\hfill \\ \multicolumn{1}{c}{}& {\left[\begin{array}{cccccccccccc}\hfill u_i^m\hfill & \hfill v_i^m\hfill & \hfill w_i^m\hfill & \hfill {\theta }_{x,i}^m\hfill & \hfill {\theta }_{y,i}^m\hfill & \hfill {\theta }_{z,i}^m\hfill & \hfill u_j^m\hfill & \hfill v_j^m\hfill & \hfill w_j^m\hfill & \hfill {\theta }_{x,j}^m\hfill & \hfill {\theta }_{y,j}^m\hfill & \hfill {\theta }_{z,j}^m\hfill \end{array}\right]}^T\hfill \end{array}$$

(2)

where $u$ , $v$ , and $w$ are the translational displacements for both ends of the beam, respectively; and $\theta$ is the angular displacement. The formula for stiffness equation can be expressed as

$$\left\{F_m\right\}=\left[k_m\right]\left\{{\delta }_m\right\}$$

(3)

where $[k_m]$ is the stiffness matrix.

A beam would disappear only when the tensile strain and bend angle exceeds permissible value. Therefore, the fracture criterion ¹⁰ is

$${\left(\frac{\epsilon }{t_{\epsilon }}\right)}^2+\frac{\mathrm{Max}(\left|{\Theta }^1\right|,\left|{\Theta }^2\right|)}{t_{\Theta }}\ge 1$$

(4)

where $\epsilon =\Delta l/l$ is axial strain of the beam ( $l$ is the initial length of the beam and $\Delta l$ is the deformation along the direction of axis), ${\Theta }^1$ and ${\Theta }^2$ are bend angles of both ends of the beam, and $t_{\epsilon }$ and $t_{\Theta }$ are threshold values for two fracture modes of extending and bending.

Two cases

Two cases were simulated to verify the validity of the method proposed in this work. The first case is to simulate a uniaxial compressive test of a rock specimen, and the second case is to simulate a blasting test of geotechnical medium with multilayers. In the two simulations, the specimens were generated by a seed expansion method (SEM). ¹⁷

The specimens used in the two simulations were generated by the SEM. ¹⁷ Seeds were first generated inside a given boundary, which include four spheres. The bases, which include three spheres, can be obtained by local Delaunay tessellation, which serve to generate new spheres by a simple algebraic calculation. This process is called as seed expansion. The distance function is employed to determine the distance between particles and boundaries. A refilling method has been developed in order to improve the packing density.

Simulation of uniaxial compressive test

A uniaxial compressive test was simulated using the method proposed in this work. The pressure load, which was generated from a compression testing machine in a uniaxial compressive test, was applied in the simulation. The specimen generated by the SEM contains 21,432 spheres. The simulated parameters are listed in Table 2. The time step is 10⁻⁶s. It took 2 h or so to implement the simulation. The simulated results are shown in Figure 3. The comparison of uniaxial compressive curves from numerical simulation and experiment is shown in Figure 4. The simulated uniaxial compressive strength is approximately 10% less than that of experiment.

Table 2.

Parameters used in simulation of uniaxial compressive test.

Parameters	Variables	Particle	Beam
Density	$\rho ({\mathrm{kg}/\mathrm{m}}^3)$	2600	–
Elastic modulus	$E(\mathrm{GPa})$	2.5	2.5
Shear modulus	$G(\mathrm{GPa})$	1.5	1.5
Poisson ratio	$v$	0.2	0.2
Normal damping	${\eta }_d$	0.3	–
Tangential damping	${\eta }_{\tau }$	0.3	–
Friction coefficient	$\mu$	0.7	–
Permissible tensile strain	$t_{\epsilon }(\%)$	–	0.5
Permissible bend angle	$t_{\theta }(\circ )$	–	1.0

Figure 3.

Configurations of uniaxial compressive test for a rock specimen (specimen size: 50 mm × 100 mm).

Figure 4.

Simulated and experimental uniaxial compressive curves.

Simulation of explosion in geotechnical media

To verify the validity of this method proposed in this work, a explosion in a geotechnical media was simulated. The specimen was also generated using the SEM. The specimen contains 405,672 spheres. For convenience, an empirical formula ¹⁸ is employed to describe the explosion pressure attenuation with time, that is

$$\Delta p_m=\frac{0.76}{\overline{R}}{\left(\frac{K}{8}\right)}^{2/3}+\frac{2.25}{{\overline{R}}^2}{\left(\frac{K}{8}\right)}^{1/3}+\frac{6.5}{{\overline{R}}^3}$$

(5)

where $\Delta p_m$ is the explosion pressure (kg/cm²), $W$ is TNT equivalence of explosives, $K$ is the bulk modulus, and $\overline{R}$ is scaled distance, that is, $\overline{R}=R/\sqrt[3]{W}$ . In this case, the value of bulk modulus $K$ is 40 GPa, the charge weight is 0.26 kg, the burial depth of the charge is 0.38 m, and the rest of the parameters used in this simulation are the same as the case one. The simulated results are shown in Figure 5(a)–(c), and Figure 5(d) is the corresponding practical explosion result. The comparison of simulated and experimental results is listed in Table 3. The results show that the simulated results approximately approach to those of practical explosion.

Table 3.

Comparison of the simulated and the explosion test results.

	Simulated results	Experimental results
Depth of explosion pit (m)	0.23	0.18
Volume of explosion pit (m3)	0.024	0.028

Figure 5.

Simulated and practical result configurations for an explosion in a geotechnical media ((a) at 0.8 ms, (b) at 2 ms, (c) at 20 ms, and (d) photograph of the practical explosion result).

Discussion and conclusion

The BPM was developed in the frame of DEM in this work. A uniaxial compressive test and explosion in technical media were simulated using the method proposed in this work. The simulated results show that (1) the fracture and fragmentation can be well modeled by the BPM and (2) the phenomena of crack generation, crack developing, fragments formation, and fragments flying can also be reflected by the simulated results. Therefore, the BPM is suitable to simulate the process of fracture and fragmentation in which large deformation and rigid motions of specimen are involved. The method proposed in this work can break through some limitations of the methods based on the continuum mechanics when simulating large deformation and rigid motions.