# Numerical Simulation Study on Flow and Fracture of Granular Materials

Authors: Mingqing Liu; Jun Liu; Guangkun Liu; Peng Shan; Zhimin Xiao; Futian Zhao
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

Flow and fracture of granular materials under external loads is a complex mechanical process, and the research on its law is still in the exploratory stage. In this paper, the flow and fracture law of granular materials is taken as the research object, and numerical algorithm compilation and program development are combined to study. Taking full advantage of the existing algorithms and developing new ones based on the existing DEM theory, a numerical simulation program for the flow and fracture of granular materials is developed. The flow and fracture process of concrete spherical granular system with diameter of 4 cm under loading rate of 70 mm/min and end of loading of 50 kN is taken as an example to verify the simulation program. At the same time, the loading experiment of the concrete spherical particle system under the same simulation conditions was also carried out. The simulation results are compared with the experimental results in three aspects: the generation location of the particle system, the relationship between the whole load and displacement, and the degree of particle breakage. The results show that the numerical simulation is in good agreement with the experimental results, which verifies the reliability of the numerical algorithm and the simulation program, and can provide support for the study of the flow and fracture process of granular materials.


---

## Body

## 1. Introduction

Granular materials are aggregates composed of a large number of discrete particles, which exist widely in nature. There is no connection between granular particles, and flow occurs under gravity or engineering loads. During the flow process, the particles in the system will break up due to the extrusion and collision between granular particles or between granular particle and boundary. The broken secondary particles will merge into the granular system and continue to flow with the system and may further break up. Granular material is ubiquitous in many engineering fields [1, 2], and its mechanical behavior is an important research topic in related fields, such as stability analysis of highly weathered rock slope in mining engineering, penetration effect of loose rock mass in protective engineering and design of shelter layer, flow and fracture of railway slag in railway engineering under the action of train pulsation, etc., so the mechanism for flow and fracture of granular materials needed to be clarified.The typical characteristics of granular materials are as follows: they are similar to solids when undisturbed, and they show fluid properties when disturbed, but they are different from solids and fluids. The concept of strain in granular materials is ambiguous, and there is no consistent description for the two states of granular materials, i.e., quasisolid and quasiliquid [3, 4]. Due to the large number of particles in granular materials and the fluidity of particle system, it is difficult to measure the contact force between particles by experimental means. With the development of computer, the numerical simulation method has become an effective way to study the mechanical behavior of granular materials.At present, the numerical simulation methods for granular materials are mainly the finite element method and discrete element method. For the finite element method, it is widely used as a numerical scheme in the continuum model [5–8]. Some scholars use the finite element method to analyze the mechanical properties of particles before they break up. When particles break up, the discrete element simulation method is used to study the process of particle breakup [9–11]. This is mainly because the simulation method based on continuum mechanics can effectively analyze the internal stress and displacement distribution characteristics of particles under external loads. However, the breaking process of particles, such as crack generation and propagation, and debris flying away from the parent body and other phenomena cannot be described.In the discrete element method, the basic element is a rigid body and cannot be broken. If a single particle in the granular system is taken as the basic element, the particle will not be broken. In fact, the stress state of the particle will change under the action of load, and there will be many tiny cracks in the particle, which will lead to the destruction of the particle. In the current numerical simulation studies, many scholars still regard particles as rigid body elements, which will lead to great differences between simulation results and actual results [12–20]. In order to overcome the defect that particles cannot be broken because they become rigid bodies, some scholars simulate the breakage of single particles [13, 14, 17, 18, 21–23] by treating single particles as aggregates of several rigid elements and adding various bonding models (such as the particle-beam model, parallel bonding model, etc.) among the elements. In two-dimensional simulation, the common element shapes are circular and polygon, while in three-dimensional simulation, spherical is the main element. However, the cluster elements consisting of rigid spheres still have voids, which are different from the real materials. Based on this, some scholars use the method of discretizing a single particle into a polyhedral element and adding beams between the elements to simulate the mechanical behavior of a particle system composed of multiple particles. Another important point is that one of the limitations of discrete element method is computational efficiency. Granular materials are aggregates of a large number of particles. If each particle is divided into several units, the number of units in the system is very large. In order to simulate the fragmentation process of granular system composed of many particles, a large amount of computing time is needed, and the performance of computer is required enormously [24]. Because of this limitation, there are few studies on discrete element simulation of multiple particle breakage, especially three-dimensional simulation, and the number of particles involved in the related studies is also relatively small.In summary, based on the existing discrete element method theory [25–30], and using the three-dimensional block element beam-particle model [31], particle failure simulation mesh generation algorithm [32], and three-dimensional block element contact detection algorithm, the simulation program of granular material flow and fracture is developed. The program can truly describe the damage evolution process of a single particle in a granular system under the action of surrounding particles or boundaries, as well as the generation and propagation of cracks and the formation of fragments and the whole process of flying away from the parent body. In order to verify the reliability of the program, the flow and fracture process of concrete spherical particle system with diameter of 4 cm under external loading was simulated by the program, and the corresponding loading experiments were carried out. By means of comparative analysis, the simulation results are in good agreement with the experimental results, which shows that the program can simulate the flow and fracture process of granular materials well, and it can provide an effective support means for numerical simulation research in this field.


## 2. Numerical Algorithms for Flow and Fracture of Granular Materials

In order to realize the development of simulation program for flow and fracture of granular materials, the numerical algorithms used are very important. The numerical algorithms needed for the simulation program are shown as follows: (1) initial configuration generation algorithm for granular system, (2) mesh generation algorithm for particle failure simulation, (3) contact detection algorithm for block elements, (4) calculation algorithm of contact force between block elements, (5) calculation method of contact force between block element and boundary, (6) calculation algorithm of motion for block elements, (7) generation and action algorithm of beam elements, (8) calculation method of global load and local load of granular particle system under external load, (9) fragmentation scale statistics algorithm, and (10) visualization algorithm. Among them, algorithm (8) is well known and algorithms (1), (4), (5), (6), and (10) can be used from the literature [28–30, 33], so this paper does not introduce them in detail. The development process and detailed information of other numerical algorithms are described below.


### 2.1. Mesh Generation Algorithm for Particle Failure Simulation

#### 2.1.1. Mesh Generation Algorithm for Single Particle

For a single particle, the algorithm in the literature [31] is used for mesh generation.


#### 2.1.2. Mesh Generation Algorithm for the Granular System

For the granular system composed of multiple particles, it will waste computing time if each particle is meshed one by one. In this paper, a particle is generated at the origin of the coordinate system. The radius of the particle can be taken as the radius of any particle in the granular system, and the particle is meshed. That is to say, the particle is discretized by the tetrahedron element, and the information of all tetrahedron elements is recorded, including the coordinates and numbers of the tetrahedron vertices. Then, the spherical center of the particle is moved to the position of the spherical center of each particle in the granular system, and all tetrahedrons are enlarged or reduced according to the ratio of particle diameter. Applying the quaternion principle [34], rotate the tetrahedrons contained in each particle in turn; the rotation axis and angle of each particle are random. Finally, update the basic information of tetrahedrons. The mesh generation algorithm for granular system is shown in Algorithm 1.Algorithm 1:Mesh generation algorithm for the particle system.
(1)
Transfer information about particle size and position in granular system obtained from initial configuration algorithm.
(2)
Generate a particle at the origin of the coordinate system. The radius of the particle can be taken as the radius of any particle in the granular system.
(3)
Start the mesh generation algorithm, discretize the particle with tetrahedron elements, and record the basic information of each tetrahedron, including the coordinates of tetrahedron vertices, which vertices make up each edge, each surface, etc.
(4)
Cycle all particles in the granular system.
(a) Translate the spherical center of the particle to the position of the spherical center of a particle in the granular system.
(b) Enlarge or reduce all tetrahedrons according to the ratio of two particle diameters.
(c) Rotate all tetrahedrons according to random rotation axis and rotation angle.
(d) Update the basic information of tetrahedrons.
(5)
End the cycle.
(6)
Generate a directory and store all the basic information of tetrahedron in the directory.


### 2.2. Contact Detection Algorithm for Block Elements

#### 2.2.1. First-Order Contact Detection Algorithm

Based on the theory of the boundary box method [35], this paper develops a first-order contact detection algorithm, which is shown in Algorithm 2.Algorithm 2:First-order contact detection algorithm.
(1)
Transfer the basic information of all tetrahedrons obtained from mesh generation algorithm of granular system.
(2)
Construct boundary box for each tetrahedron.
(3)
Project the boundary box onto three coordinate axes of granular system.
(4)
Generate three catalogues and store the sequence of starting points and ending points of boundary box projecting on three coordinate axes in the catalogue.
(5)
Cycle all tetrahedrons (for (i = 0; i < num; i++) (num is the number of tetrahedrons)).
(a) Find the starting point and ending point of the projection of the boundary box of the tetrahedroni on the three coordinate axes from the storage catalogue.
(b) Cycle all remaining tetrahedrons after removing the tetrahedroni (for (j = 0; j < num-1; i++)).
① Generate a directory.
② Judge whether the projection on three coordinate axes of the boundary box of tetrahedronj and tetrahedron i are overlapping; if overlap, then the tetrahedron j is considered as the neighbor unit of the tetrahedron i, and the unit is stored in the catalogue.
(c) End the cycle.
(6)
End the cycle.
(7)
Get the directory of all tetrahedron neighborhood units.


#### 2.2.2. Second-Order Contact Detection Algorithm

In this paper, the fast common-plane method [36] is selected to determine whether the blocks are in contact with each other and obtain detailed contact information. On the basis of this theory, a second-order contact detection algorithm for block elements is developed, as shown in Algorithm 3.Algorithm 3:Second-order contact detection algorithm.
(1)
Transfer geometric information of two tetrahedron neighbor units A and B that need to determine whether there is contact between them.
(2)
Define a set of vecN of three-dimensional vectors to store the unit normal vector of the common surface in each iteration step.
(3)
Connect the centroid of A and B and make the vertical bisector of the connecting line. With the vertical bisector as the common surface, the nearest verticesvA and vB. of A and B are obtained.
(4)
Calculate the direction vectors of the connecting line ofvA and vB; the direction vector is stored in vecN as the unit normal vector of the common surface.
(5)
The faces sharing verticesvA and vB are regarded as potential common surfaces and their unit normal vectors are stored in vecN, respectively.
(6)
Make vector products for the direction vectors of the edges sharingvA vertex and vB vertex, respectively, unitized the vector product, and it is also stored in vecN as a normal vector of potential common surface.
(7)
Take an edge that contains vertexvA, calculate its direction vector, and then calculate the direction vector of the line between vA and vB; calculate their vector product, then take it to vector product with the direction vector of the edge passing through vA vertex, and store the vector product in vecN as the normal vector of the potential common surface. Repeat the above operations for all edges passing through vertex vA, and then repeat the above operations for all edges passing through vertex vB.
(8)
For the normal vectors of each potential common surface in the set vecN, calculate the minimum distance between the two tetrahedrons, marking the nearest vertex asvA′ and vB′; if vA′=vA and vB′=vB, then the normal vector of the common surface is the final normal vector of the common surface and terminates the calculation; otherwise, return to 4 and continue the calculation.


### 2.3. Generation and Action Algorithm of Beam Elements

The algorithm of beam element generation and action adopted in this paper are shown in Algorithms4 and 5, respectively.Algorithm 4:Beam element generation algorithm.
(1)
Obtain the detailed contact information of tetrahedron according to the algorithm in Section2.2.
(2)
Cycle for all tetrahedrons (for (i = 0; i < num; i++) (num is the number of tetrahedrons)).
(a) Cycle all tetrahedrons in contact with tetrahedronj (for (j = 0; j < neighborNum; j++) (neighborNum is the number of tetrahedrons in contact with tetrahedron i)).
If there is no beam between tetrahedroni and tetrahedron j.
① Get the centroid of tetrahedroni and neighbor unit j.
② The starting points and ending points of the beam are set to be the mass centers of two tetrahedrons, respectively, and the length of the beam is the distance between the centers of mass.
③ Set the cross-sectional area of the beam as the contact area.
④ Girder is endowed with mechanical properties parameters.
(b) End the cycle.
(3)
End the cycle.
(4)
Generate the catalogue of interblock beams.Algorithm 5:Beam element action algorithm.
(1)
The list of beams between all block elements is obtained from the algorithm in Section2.2.
(2)
Cycle all the beams in the catalogue (for (i = 0; i < beamNum; i++) (beamNum is the number of beam)).
(a) Obtain information about two blocks connected by beami.
(b) Obtain the displacement and rotation angles of two blocks.
(c) Calculate the relative displacement and rotation angle at both ends of the beam.
(d) Deformation of the beam is obtained by comparing the position of the beam with that of the previous step-length beam.
(e) Calculation of forces, moments, and torques on beams based on the matrix displacement method.
(f) Judge by failure criteria whether the beam disappears. If it disappears, then delete the beam from the catalogue.
(3)
End the cycle.


### 2.4. Fragmentation Scale Statistics Algorithm

The scale distribution of debris is an important criterion for measuring particle breakage after the particle is broken. In order to study the factors and the degree of particle breakage, a statistical algorithm of debris size was developed to calculate the scale distribution of debris in numerical simulation.In the fragment size statistics algorithm, the first thing is to identify which tetrahedron the fragments is composed of, then obtain the position and size information of all the tetrahedrons that make up the fragment, determine the shape and size of the fragments according to these information, and at last, judge whether the fragments can pass through the sieve hole or not. After screening, the volume of the fragments of each sieve is calculated in turn; dividing by the sum of the volume of all particles in the granular system, the percentage of the mass of the fragments smaller than a certain particle size to the total mass can be obtained, and the scale distribution of the fragments can be obtained. The fragment screening algorithm is shown in Algorithm6.Algorithm 6:Fragmentation sieving algorithm.
(1)
Transfer information of all tetrahedrons and beams in granular system after loading.
(2)
Generate fragmentation directories. Each fragment in the catalogue consists of the number of tetrahedrons it contains, and the initial is empty.
(3)
Cycle all tetrahedrons (for (i = 0; i < num; i++) (num is the number of tetrahedrons)).
(a) Determine whether the generated fragments contain tetrahedroni.
① If included.
Continue with the next tetrahedron cycle.
② Otherwise, judge whether there is a beam between tetrahedroni and other tetrahedrons.
(i) If existed.
(1) Record the tetrahedron number and continue to determine whether the tetrahedron is associated with other tetrahedrons (excluding tetrahedroni), until there is no new tetrahedron.
(2) Generate fragments; the fragments consist of all tetrahedron numbers found in this search.
(ii) Otherwise.
Record the number of the tetrahedron and mark the tetrahedron as a fragment.
(b) Update the fragmentation directory.
(4)
End the cycle.
(5)
Cycle all sieves with all aperture (for (j = 0; j < screenNum; j++) (screenNum is the number of sieve)).
(a) Cycle all debris (for (k = 0; k < fragNum; k++) (fragNum is the number of fragments)). According to the position and size information of all tetrahedrons contained in fragment k, judge whether fragment k can pass through sieve and record fragmentation information if passed.
(b) Calculate the sum of the volumes of all fragments passing through sievej.
(c) Calculate the percentage of the mass of fragments smaller than the sievej aperture to the total mass.
(d) End the cycle.
(6)
End the cycle.


## 2.1. Mesh Generation Algorithm for Particle Failure Simulation

### 2.1.1. Mesh Generation Algorithm for Single Particle

For a single particle, the algorithm in the literature [31] is used for mesh generation.


### 2.1.2. Mesh Generation Algorithm for the Granular System

For the granular system composed of multiple particles, it will waste computing time if each particle is meshed one by one. In this paper, a particle is generated at the origin of the coordinate system. The radius of the particle can be taken as the radius of any particle in the granular system, and the particle is meshed. That is to say, the particle is discretized by the tetrahedron element, and the information of all tetrahedron elements is recorded, including the coordinates and numbers of the tetrahedron vertices. Then, the spherical center of the particle is moved to the position of the spherical center of each particle in the granular system, and all tetrahedrons are enlarged or reduced according to the ratio of particle diameter. Applying the quaternion principle [34], rotate the tetrahedrons contained in each particle in turn; the rotation axis and angle of each particle are random. Finally, update the basic information of tetrahedrons. The mesh generation algorithm for granular system is shown in Algorithm 1.Algorithm 1:Mesh generation algorithm for the particle system.
(1)
Transfer information about particle size and position in granular system obtained from initial configuration algorithm.
(2)
Generate a particle at the origin of the coordinate system. The radius of the particle can be taken as the radius of any particle in the granular system.
(3)
Start the mesh generation algorithm, discretize the particle with tetrahedron elements, and record the basic information of each tetrahedron, including the coordinates of tetrahedron vertices, which vertices make up each edge, each surface, etc.
(4)
Cycle all particles in the granular system.
(a) Translate the spherical center of the particle to the position of the spherical center of a particle in the granular system.
(b) Enlarge or reduce all tetrahedrons according to the ratio of two particle diameters.
(c) Rotate all tetrahedrons according to random rotation axis and rotation angle.
(d) Update the basic information of tetrahedrons.
(5)
End the cycle.
(6)
Generate a directory and store all the basic information of tetrahedron in the directory.


## 2.1.1. Mesh Generation Algorithm for Single Particle

For a single particle, the algorithm in the literature [31] is used for mesh generation.


## 2.1.2. Mesh Generation Algorithm for the Granular System

For the granular system composed of multiple particles, it will waste computing time if each particle is meshed one by one. In this paper, a particle is generated at the origin of the coordinate system. The radius of the particle can be taken as the radius of any particle in the granular system, and the particle is meshed. That is to say, the particle is discretized by the tetrahedron element, and the information of all tetrahedron elements is recorded, including the coordinates and numbers of the tetrahedron vertices. Then, the spherical center of the particle is moved to the position of the spherical center of each particle in the granular system, and all tetrahedrons are enlarged or reduced according to the ratio of particle diameter. Applying the quaternion principle [34], rotate the tetrahedrons contained in each particle in turn; the rotation axis and angle of each particle are random. Finally, update the basic information of tetrahedrons. The mesh generation algorithm for granular system is shown in Algorithm 1.Algorithm 1:Mesh generation algorithm for the particle system.
(1)
Transfer information about particle size and position in granular system obtained from initial configuration algorithm.
(2)
Generate a particle at the origin of the coordinate system. The radius of the particle can be taken as the radius of any particle in the granular system.
(3)
Start the mesh generation algorithm, discretize the particle with tetrahedron elements, and record the basic information of each tetrahedron, including the coordinates of tetrahedron vertices, which vertices make up each edge, each surface, etc.
(4)
Cycle all particles in the granular system.
(a) Translate the spherical center of the particle to the position of the spherical center of a particle in the granular system.
(b) Enlarge or reduce all tetrahedrons according to the ratio of two particle diameters.
(c) Rotate all tetrahedrons according to random rotation axis and rotation angle.
(d) Update the basic information of tetrahedrons.
(5)
End the cycle.
(6)
Generate a directory and store all the basic information of tetrahedron in the directory.


## 2.2. Contact Detection Algorithm for Block Elements

### 2.2.1. First-Order Contact Detection Algorithm

Based on the theory of the boundary box method [35], this paper develops a first-order contact detection algorithm, which is shown in Algorithm 2.Algorithm 2:First-order contact detection algorithm.
(1)
Transfer the basic information of all tetrahedrons obtained from mesh generation algorithm of granular system.
(2)
Construct boundary box for each tetrahedron.
(3)
Project the boundary box onto three coordinate axes of granular system.
(4)
Generate three catalogues and store the sequence of starting points and ending points of boundary box projecting on three coordinate axes in the catalogue.
(5)
Cycle all tetrahedrons (for (i = 0; i < num; i++) (num is the number of tetrahedrons)).
(a) Find the starting point and ending point of the projection of the boundary box of the tetrahedroni on the three coordinate axes from the storage catalogue.
(b) Cycle all remaining tetrahedrons after removing the tetrahedroni (for (j = 0; j < num-1; i++)).
① Generate a directory.
② Judge whether the projection on three coordinate axes of the boundary box of tetrahedronj and tetrahedron i are overlapping; if overlap, then the tetrahedron j is considered as the neighbor unit of the tetrahedron i, and the unit is stored in the catalogue.
(c) End the cycle.
(6)
End the cycle.
(7)
Get the directory of all tetrahedron neighborhood units.


### 2.2.2. Second-Order Contact Detection Algorithm

In this paper, the fast common-plane method [36] is selected to determine whether the blocks are in contact with each other and obtain detailed contact information. On the basis of this theory, a second-order contact detection algorithm for block elements is developed, as shown in Algorithm 3.Algorithm 3:Second-order contact detection algorithm.
(1)
Transfer geometric information of two tetrahedron neighbor units A and B that need to determine whether there is contact between them.
(2)
Define a set of vecN of three-dimensional vectors to store the unit normal vector of the common surface in each iteration step.
(3)
Connect the centroid of A and B and make the vertical bisector of the connecting line. With the vertical bisector as the common surface, the nearest verticesvA and vB. of A and B are obtained.
(4)
Calculate the direction vectors of the connecting line ofvA and vB; the direction vector is stored in vecN as the unit normal vector of the common surface.
(5)
The faces sharing verticesvA and vB are regarded as potential common surfaces and their unit normal vectors are stored in vecN, respectively.
(6)
Make vector products for the direction vectors of the edges sharingvA vertex and vB vertex, respectively, unitized the vector product, and it is also stored in vecN as a normal vector of potential common surface.
(7)
Take an edge that contains vertexvA, calculate its direction vector, and then calculate the direction vector of the line between vA and vB; calculate their vector product, then take it to vector product with the direction vector of the edge passing through vA vertex, and store the vector product in vecN as the normal vector of the potential common surface. Repeat the above operations for all edges passing through vertex vA, and then repeat the above operations for all edges passing through vertex vB.
(8)
For the normal vectors of each potential common surface in the set vecN, calculate the minimum distance between the two tetrahedrons, marking the nearest vertex asvA′ and vB′; if vA′=vA and vB′=vB, then the normal vector of the common surface is the final normal vector of the common surface and terminates the calculation; otherwise, return to 4 and continue the calculation.


## 2.2.1. First-Order Contact Detection Algorithm

Based on the theory of the boundary box method [35], this paper develops a first-order contact detection algorithm, which is shown in Algorithm 2.Algorithm 2:First-order contact detection algorithm.
(1)
Transfer the basic information of all tetrahedrons obtained from mesh generation algorithm of granular system.
(2)
Construct boundary box for each tetrahedron.
(3)
Project the boundary box onto three coordinate axes of granular system.
(4)
Generate three catalogues and store the sequence of starting points and ending points of boundary box projecting on three coordinate axes in the catalogue.
(5)
Cycle all tetrahedrons (for (i = 0; i < num; i++) (num is the number of tetrahedrons)).
(a) Find the starting point and ending point of the projection of the boundary box of the tetrahedroni on the three coordinate axes from the storage catalogue.
(b) Cycle all remaining tetrahedrons after removing the tetrahedroni (for (j = 0; j < num-1; i++)).
① Generate a directory.
② Judge whether the projection on three coordinate axes of the boundary box of tetrahedronj and tetrahedron i are overlapping; if overlap, then the tetrahedron j is considered as the neighbor unit of the tetrahedron i, and the unit is stored in the catalogue.
(c) End the cycle.
(6)
End the cycle.
(7)
Get the directory of all tetrahedron neighborhood units.


## 2.2.2. Second-Order Contact Detection Algorithm

In this paper, the fast common-plane method [36] is selected to determine whether the blocks are in contact with each other and obtain detailed contact information. On the basis of this theory, a second-order contact detection algorithm for block elements is developed, as shown in Algorithm 3.Algorithm 3:Second-order contact detection algorithm.
(1)
Transfer geometric information of two tetrahedron neighbor units A and B that need to determine whether there is contact between them.
(2)
Define a set of vecN of three-dimensional vectors to store the unit normal vector of the common surface in each iteration step.
(3)
Connect the centroid of A and B and make the vertical bisector of the connecting line. With the vertical bisector as the common surface, the nearest verticesvA and vB. of A and B are obtained.
(4)
Calculate the direction vectors of the connecting line ofvA and vB; the direction vector is stored in vecN as the unit normal vector of the common surface.
(5)
The faces sharing verticesvA and vB are regarded as potential common surfaces and their unit normal vectors are stored in vecN, respectively.
(6)
Make vector products for the direction vectors of the edges sharingvA vertex and vB vertex, respectively, unitized the vector product, and it is also stored in vecN as a normal vector of potential common surface.
(7)
Take an edge that contains vertexvA, calculate its direction vector, and then calculate the direction vector of the line between vA and vB; calculate their vector product, then take it to vector product with the direction vector of the edge passing through vA vertex, and store the vector product in vecN as the normal vector of the potential common surface. Repeat the above operations for all edges passing through vertex vA, and then repeat the above operations for all edges passing through vertex vB.
(8)
For the normal vectors of each potential common surface in the set vecN, calculate the minimum distance between the two tetrahedrons, marking the nearest vertex asvA′ and vB′; if vA′=vA and vB′=vB, then the normal vector of the common surface is the final normal vector of the common surface and terminates the calculation; otherwise, return to 4 and continue the calculation.


## 2.3. Generation and Action Algorithm of Beam Elements

The algorithm of beam element generation and action adopted in this paper are shown in Algorithms4 and 5, respectively.Algorithm 4:Beam element generation algorithm.
(1)
Obtain the detailed contact information of tetrahedron according to the algorithm in Section2.2.
(2)
Cycle for all tetrahedrons (for (i = 0; i < num; i++) (num is the number of tetrahedrons)).
(a) Cycle all tetrahedrons in contact with tetrahedronj (for (j = 0; j < neighborNum; j++) (neighborNum is the number of tetrahedrons in contact with tetrahedron i)).
If there is no beam between tetrahedroni and tetrahedron j.
① Get the centroid of tetrahedroni and neighbor unit j.
② The starting points and ending points of the beam are set to be the mass centers of two tetrahedrons, respectively, and the length of the beam is the distance between the centers of mass.
③ Set the cross-sectional area of the beam as the contact area.
④ Girder is endowed with mechanical properties parameters.
(b) End the cycle.
(3)
End the cycle.
(4)
Generate the catalogue of interblock beams.Algorithm 5:Beam element action algorithm.
(1)
The list of beams between all block elements is obtained from the algorithm in Section2.2.
(2)
Cycle all the beams in the catalogue (for (i = 0; i < beamNum; i++) (beamNum is the number of beam)).
(a) Obtain information about two blocks connected by beami.
(b) Obtain the displacement and rotation angles of two blocks.
(c) Calculate the relative displacement and rotation angle at both ends of the beam.
(d) Deformation of the beam is obtained by comparing the position of the beam with that of the previous step-length beam.
(e) Calculation of forces, moments, and torques on beams based on the matrix displacement method.
(f) Judge by failure criteria whether the beam disappears. If it disappears, then delete the beam from the catalogue.
(3)
End the cycle.


## 2.4. Fragmentation Scale Statistics Algorithm

The scale distribution of debris is an important criterion for measuring particle breakage after the particle is broken. In order to study the factors and the degree of particle breakage, a statistical algorithm of debris size was developed to calculate the scale distribution of debris in numerical simulation.In the fragment size statistics algorithm, the first thing is to identify which tetrahedron the fragments is composed of, then obtain the position and size information of all the tetrahedrons that make up the fragment, determine the shape and size of the fragments according to these information, and at last, judge whether the fragments can pass through the sieve hole or not. After screening, the volume of the fragments of each sieve is calculated in turn; dividing by the sum of the volume of all particles in the granular system, the percentage of the mass of the fragments smaller than a certain particle size to the total mass can be obtained, and the scale distribution of the fragments can be obtained. The fragment screening algorithm is shown in Algorithm6.Algorithm 6:Fragmentation sieving algorithm.
(1)
Transfer information of all tetrahedrons and beams in granular system after loading.
(2)
Generate fragmentation directories. Each fragment in the catalogue consists of the number of tetrahedrons it contains, and the initial is empty.
(3)
Cycle all tetrahedrons (for (i = 0; i < num; i++) (num is the number of tetrahedrons)).
(a) Determine whether the generated fragments contain tetrahedroni.
① If included.
Continue with the next tetrahedron cycle.
② Otherwise, judge whether there is a beam between tetrahedroni and other tetrahedrons.
(i) If existed.
(1) Record the tetrahedron number and continue to determine whether the tetrahedron is associated with other tetrahedrons (excluding tetrahedroni), until there is no new tetrahedron.
(2) Generate fragments; the fragments consist of all tetrahedron numbers found in this search.
(ii) Otherwise.
Record the number of the tetrahedron and mark the tetrahedron as a fragment.
(b) Update the fragmentation directory.
(4)
End the cycle.
(5)
Cycle all sieves with all aperture (for (j = 0; j < screenNum; j++) (screenNum is the number of sieve)).
(a) Cycle all debris (for (k = 0; k < fragNum; k++) (fragNum is the number of fragments)). According to the position and size information of all tetrahedrons contained in fragment k, judge whether fragment k can pass through sieve and record fragmentation information if passed.
(b) Calculate the sum of the volumes of all fragments passing through sievej.
(c) Calculate the percentage of the mass of fragments smaller than the sievej aperture to the total mass.
(d) End the cycle.
(6)
End the cycle.


## 3. Programming for Flow and Fracture of Granular Materials

This paper chooses Visual Studio .NET 2003 as the development platform and makes full use of STL resource base and C++ language, and based on the numerical algorithms of the second part, the numerical simulation program of granular materials flow and fracture is developed, and the visualization of simulation results is realized by using OpenGL technology.The developed simulation program adopts modular design. And the advantage of using modular programming is that complex program can be decomposed into several simple programs. The program is readable and easy to modify, improves efficiency, and saves time for program development. The whole program consists of three modules: preprocessing module, main calculation module, and visualization module, as shown in Figure1.Figure 1
Function module division of flow and fracture simulation program for granular materials.


### 3.1. Preprocessing Module

The main functions of the preprocessing module are shown as follows: (1) read in the geometric characteristics of container for granular materials needed for the simulation, such as, the geometric characteristics and mechanical parameters of the pressure plate, the coordinates of the measuring points, and the mechanical parameters of granular particles, (2) generate granular system based on the initial configuration generation algorithm, and (3) according to the mesh generation algorithm of granular system, mesh all granular particles. That is to say, the particle is discretized by tetrahedron elements. And generate data files of tetrahedron units to provide initial data for subsequent calculations. The flow chart of the preprocessing module is shown in Figure2.Figure 2
Flow chart for preprocessing module program.


### 3.2. Main Calculation Module

The main calculation module is the core content of the developed simulation program. Its main functions include the following: (1) Judge whether there is contact between blocks based on contact detection algorithm of block elements and calculate detailed contact information such as contact point, contact depth, contact area, and contact normal direction. (2) Calculate interblock contact force based on the calculation method of interblock contact force. (3) Calculate contact force between block element and boundary based on contact force calculation algorithm between block element and boundary. (4) Update the acceleration, velocity, position, and rotation angle of the block based on the block element motion algorithm. (5) Generate beams between contacted block elements based on generation action algorithm of beam elements. Calculate the mechanical responses of beams such as end force and deformation of beams based on relative displacement and rotation angle between block elements and judge whether the beam disappears or not based on the strength criterion of the beam. (6) Output the basic information of all tetrahedrons and the beam in each calculation step in the file format to provide data for the visualization of the simulation results. The program flow chart of the main calculation module is shown in Figure3.Figure 3
Flow chart of the visualization module program.


### 3.3. Visualization Module

The visualization module makes full use of the efficient STL (Standard Template Library, Standard Template Library) resource library, and the visualization program is compiled with OpenGL graphics library. Visualization module has graphics output function, and drawing commands are given in the form of data files (parameter tables). When the program runs, it draws all kinds of information recorded in the process of structure destruction one by one according to different synchronous data files and stores the pictures on disk in order. The flow chart of the visualization module program is shown in Figure4.Figure 4
Flow chart for visualization module program.


### 3.4. Verification for Program

In order to verify the reliability of the numerical algorithm and the simulation program developed in this paper, a numerical example of single particle falling naturally from a height of 20 cm and impacting the ground at an initial velocity of 20 m/s is simulated. And experiments of single concrete sphere impacting the ground under the same conditions were carried out. The simulation results and experimental results are shown in Figures5–9.Figure 5
Simulation results of single particle falling naturally to the ground: (a) tetrahedron elements; (b) beam grid.
(a)
(b)Figure 6
Experimental results of single concrete sphere falling to the ground from a height of 20 cm.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)Figure 7
Simulation results of single particle falling to the ground at an initial velocity of 20 m/s.Figure 8
Beam failure process of single particle impacting ground at an initial velocity of 20 m/s.
(a)
(b)
(c)
(d)Figure 9
Experimental results of single concrete sphere falling to the ground at an initial velocity of 20 m/s.
(a)
(b)
(c)
(d)
(e)
(f)It can be seen from Figure5(a) that the particle falls naturally from a height of 20 cm. When the particle collides with the ground (the contact point is the black point in the figure), the particle bounces slightly, and the tetrahedral elements in the particle do not deform obviously. From Figure 5(b), it can be seen that the beams do not disappear after the particle collides with the ground; it indicates that the particle has not been destroyed.It can be seen from Figure6 that the concrete sphere falls freely from a height of 20 cm (the red dot is the first contact point between the sphere and the ground) and rebounds after contacting with the ground. During the contact, the sphere does not have obvious deformation or breakage, which is in good agreement with the simulation results of a single particle.As can be seen from Figures7 and 8, when the particle falls to the ground at an initial velocity of 20 m/s, a large number of beams at the bottom of the particle disappear, forming a lot of small fragments. At the same time, cracks occur and propagate in the upper part of the particle, which eventually break into several larger fragments. The simulation results are basically consistent with the actual situation.As can be seen from Figure9, single concrete sphere drops and strikes the ground at an initial velocity of 20 m/s. At the moment of contact between the sphere and the ground, the bottom of the sphere begins to break (Figure 9(b)), and then fragments are generated and fly out (Figures 9(c)–9(e)). The final breakage of the pellets is shown in Figure 9(f), in which the larger size of the pellets is the majority. The correctness of the simulation example is further verified.Through the analysis of the above simulation results and the experimental results, it can be seen that the simulation results are basically in agreement with the experimental results, which verifies the reliability of the algorithm and the program.


## 3.1. Preprocessing Module

The main functions of the preprocessing module are shown as follows: (1) read in the geometric characteristics of container for granular materials needed for the simulation, such as, the geometric characteristics and mechanical parameters of the pressure plate, the coordinates of the measuring points, and the mechanical parameters of granular particles, (2) generate granular system based on the initial configuration generation algorithm, and (3) according to the mesh generation algorithm of granular system, mesh all granular particles. That is to say, the particle is discretized by tetrahedron elements. And generate data files of tetrahedron units to provide initial data for subsequent calculations. The flow chart of the preprocessing module is shown in Figure2.Figure 2
Flow chart for preprocessing module program.


## 3.2. Main Calculation Module

The main calculation module is the core content of the developed simulation program. Its main functions include the following: (1) Judge whether there is contact between blocks based on contact detection algorithm of block elements and calculate detailed contact information such as contact point, contact depth, contact area, and contact normal direction. (2) Calculate interblock contact force based on the calculation method of interblock contact force. (3) Calculate contact force between block element and boundary based on contact force calculation algorithm between block element and boundary. (4) Update the acceleration, velocity, position, and rotation angle of the block based on the block element motion algorithm. (5) Generate beams between contacted block elements based on generation action algorithm of beam elements. Calculate the mechanical responses of beams such as end force and deformation of beams based on relative displacement and rotation angle between block elements and judge whether the beam disappears or not based on the strength criterion of the beam. (6) Output the basic information of all tetrahedrons and the beam in each calculation step in the file format to provide data for the visualization of the simulation results. The program flow chart of the main calculation module is shown in Figure3.Figure 3
Flow chart of the visualization module program.


## 3.3. Visualization Module

The visualization module makes full use of the efficient STL (Standard Template Library, Standard Template Library) resource library, and the visualization program is compiled with OpenGL graphics library. Visualization module has graphics output function, and drawing commands are given in the form of data files (parameter tables). When the program runs, it draws all kinds of information recorded in the process of structure destruction one by one according to different synchronous data files and stores the pictures on disk in order. The flow chart of the visualization module program is shown in Figure4.Figure 4
Flow chart for visualization module program.


## 3.4. Verification for Program

In order to verify the reliability of the numerical algorithm and the simulation program developed in this paper, a numerical example of single particle falling naturally from a height of 20 cm and impacting the ground at an initial velocity of 20 m/s is simulated. And experiments of single concrete sphere impacting the ground under the same conditions were carried out. The simulation results and experimental results are shown in Figures5–9.Figure 5
Simulation results of single particle falling naturally to the ground: (a) tetrahedron elements; (b) beam grid.
(a)
(b)Figure 6
Experimental results of single concrete sphere falling to the ground from a height of 20 cm.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)Figure 7
Simulation results of single particle falling to the ground at an initial velocity of 20 m/s.Figure 8
Beam failure process of single particle impacting ground at an initial velocity of 20 m/s.
(a)
(b)
(c)
(d)Figure 9
Experimental results of single concrete sphere falling to the ground at an initial velocity of 20 m/s.
(a)
(b)
(c)
(d)
(e)
(f)It can be seen from Figure5(a) that the particle falls naturally from a height of 20 cm. When the particle collides with the ground (the contact point is the black point in the figure), the particle bounces slightly, and the tetrahedral elements in the particle do not deform obviously. From Figure 5(b), it can be seen that the beams do not disappear after the particle collides with the ground; it indicates that the particle has not been destroyed.It can be seen from Figure6 that the concrete sphere falls freely from a height of 20 cm (the red dot is the first contact point between the sphere and the ground) and rebounds after contacting with the ground. During the contact, the sphere does not have obvious deformation or breakage, which is in good agreement with the simulation results of a single particle.As can be seen from Figures7 and 8, when the particle falls to the ground at an initial velocity of 20 m/s, a large number of beams at the bottom of the particle disappear, forming a lot of small fragments. At the same time, cracks occur and propagate in the upper part of the particle, which eventually break into several larger fragments. The simulation results are basically consistent with the actual situation.As can be seen from Figure9, single concrete sphere drops and strikes the ground at an initial velocity of 20 m/s. At the moment of contact between the sphere and the ground, the bottom of the sphere begins to break (Figure 9(b)), and then fragments are generated and fly out (Figures 9(c)–9(e)). The final breakage of the pellets is shown in Figure 9(f), in which the larger size of the pellets is the majority. The correctness of the simulation example is further verified.Through the analysis of the above simulation results and the experimental results, it can be seen that the simulation results are basically in agreement with the experimental results, which verifies the reliability of the algorithm and the program.


## 4. Numerical Simulation and Experiment for Flow and Fracture of Granular Materials

In order to further verify the validity of the developed numerical algorithm and program for simulating the flow and fracture of granular materials, the flow and fracture process of concrete spherical particle system with diameter of 4 cm was simulated under conditions that the loading rate is 70 mm/min and the loading will completed when the loading value is 50 kN. The numerical simulation schematic diagram is shown in Figure10. The flow and fracture experiments of concrete spherical particle system were carried out under conditions the same to the numerical simulation. And the numerical simulation results and experimental results were analyzed through comparison. The main contents of the comparison were the generation of the particle system, the relationship between the whole load and displacement in the loading process, and the size distribution of the fragments.Figure 10
Diagram of numerical simulation.


### 4.1. Numerical Simulation for Flow and Fracture of Granular Materials

#### 4.1.1. Numerical Sample Generation and Mesh Generation

The method for generating numerical samples of granular materials is as follows: the bottom particles are generated according to particle size and location information, and the particles are added one by one in a certain order. Each time the particles are added, the particles will fall naturally under the acceleration of gravity. Then, the contact force between particles and the contact force between particles and boundary are calculated according to the contact information. Particles move under the action of unbalanced force. The acceleration, velocity, and displacement of particles are calculated according to the resultant force of particles. The positions of particles are adjusted until the particle system reaches equilibrium. The generated numerical sample is shown in Figures11(d)–11(f), and the particle meshing is shown in Figure 12.Figure 11
Comparisons of particle positions recorded in experiment of concrete spherical particle system with 4 cm diameter with those generated in numerical samples: (a–c) particle positions recorded in experiment; (d–f) particle positions in numerical samples.
(a)
(b)
(c)
(d)
(e)
(f)Figure 12
Particle mesh generation: (a) single particle; (b) multiple particles.
(a)
(b)


#### 4.1.2. Process for Flow and Fracture of the Particle System

The concrete spherical particle system with diameter of 4 cm is loaded to 50 kN. The process of particle flow and fracture is shown in Figure13, in which the particle system has three layers from up to down like the situation in Figure 10, and there are 4 particles in each layer.Figure 13
Process for flow and fracture of concrete spherical particle system with 4 cm diameter.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)


#### 4.1.3. Beam Failure Process

The disappearance process of beams between tetrahedral elements is shown in Figure14 when the concrete spherical particle system with diameter of 4 cm is loaded to 50 kN.Figure 14
Disappearance process of beams between tetrahedral elements of concrete spherical particle system with 4 cm diameter.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)From Figure14, it can be seen that when the pressure plate contacts the top particles, the top particles break up and the beams between tetrahedral elements disappear. As the load increases, the number of missing beams increases. When the loading is completed, most of the beams disappear and many separate tetrahedral elements are formed, which indicates that a large amount of powder debris is produced after the particles are broken. The number of beams disappeared in the top and bottom particles were obviously more than that in the middle layer. Most of the debris in the top and bottom layer contained only two or three beams, and the beam network composed of multiple beams could be clearly seen in the middle layer, which indicated that the degree of broken particles in the top and bottom layer was more serious than that in the middle layer.


### 4.2. Experimental Verification

In order to verify the validity of the simulation program more intuitively, experiments on flow and fracture of concrete spherical particle system are carried out under the conditions the same to the simulation.The layering position of the experimental particles and the layering and crushing of the particles after loading are shown in Figure15, and the layering sieving curve and the whole sieving curve of the particle system after loading are shown in Figure 16. The whole load-displacement curve and the local load-displacement curve of the granular system during loading are shown in Figure 17.Figure 15
Particle placement before loading and the particle stratification and fragmentation after loading in the experiment: (a–c) particle placement; (d–f) surface fragmentation of each layer; (g–i) fragmentation of each layer after taking out the layers.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)Figure 16
Layered sieving curve and overall sieving curve of particle system after experimental loading.Figure 17
Whole load-displacement relationship curve and local load-displacement relationship curve of experimental granular system: (a) whole load-displacement curve; (b) local load-displacement curve.
(a)
(b)From Figures15 and 16, it can be seen that the surface particles of 4 cm concrete spherical particle system are broken most seriously under 50 kN load. Most of the fragments are in powder shape. The mass of fragments with a size less than 1 mm accounts for 28% of the total mass of the layer, the fragments with a size less than 10 mm account for 59%, and the maximum particle size of the fragments is less than 25 mm. The fragmentation degree of the lower particles is relatively light, and larger fragments can also be seen. It is found that the mass of fragments with a size less than 1 mm accounts for 24% of the total mass of the layer, the fragments with a size less than 10 mm accounts for 41%, and the maximum particle size of the fragments exceeds 45 mm. The fragmentation degree of the intermediate layer is the lightest, and the mass of fragments with a size less than 10 mm accounts for only 5% of the total mass of the layer. The maximum particle size of the fragments exceeds 45 mm. Most of the particles remain in their original shape, some of them only produce cracks or have smaller fragments peeling off the particles, and some of them break up. The average size of the fragments is much larger than that of the upper and lower layers.In the experiment, the particle system is loaded according to the displacement of the pressure plate. The loading rate is 70 mm/min, and the loading value of the granular system can be read by pushing machine during the loading process. When the load value reaches 50 kN, the loading ends. From Figure17(a), it can be seen that the whole load-displacement relationship curve presents a zigzag shape, to rise gradually in a fluctuating manner, and the curve fluctuates greatly in the initial stage. When the displacement of the pressure plate reaches 6 cm, the fluctuation range of the curve decreases obviously, and the increasing rate of load also increases significantly. When the load value is 50 kN, the displacement of the particle system is 6.8 cm. As can be seen from Figure 17(b), the local load-displacement curves obtained from the two monitoring points are quite different in both trend and magnitude; the reason for this difference is mainly related to the gradation and arrangement of the particle system. The particle system is arranged in a random manner. Under the action of external loads, the force chains with different strength will be formed inside the particle system so that the forces on the bottom of the particle system are not uniform, and the loads are not distributed equally according to the number of contact points, which is significantly different from the forces on continuous media.


### 4.3. Comparison between the Simulation Results and the Experimental Results

The comparison between the experimental results and the simulation results includes the following three aspects: the position of particle system formation, the relationship between whole load and displacement, and the scale distribution of fragments.


#### 4.3.1. Particle System Formation

The experimental photographs of particle positions in a 4 cm diameter concrete spherical particle system are compared with the generated numerical samples as shown in Figure11.From Figure11, it can be seen that the position of particles in the numerical samples of the particle system generated by the preprocessing program is basically the same as that recorded in the experiment. The particles generated by the program are standard spheres, while the concrete spheres in the experiment are manufactured by the die. The regularity of the spheres is different from that of the particles generated by the numerical simulation, but the error is not large. The generated particle system can meet the requirements of simulating granular materials’ flow and fracture under the same conditions.


#### 4.3.2. Relationship between Whole Load and Displacement

During the experiment, the particle system was loaded by moving the position of the pressure plate, and the descending speed of the pressure plate was 70 mm/min, which was consistent with the numerical simulation loading rate. In the process of loading the particle system, the contact forces between the tetrahedral elements and the pressure plate are calculated separately. The total contact forces are summed up as the resultant forces of the pressure plate. Finally, the component of the resultant forces perpendicular to the direction of the pressure plate is calculated, and the whole load of the particle system can be obtained. Based on the displacement of the pressure plate and the whole load of the particle system, the relationship curve between the whole load and the displacement of the particle is drawn and compared with the experimental results. The comparison of the relationship curve between the whole load and the displacement of concrete spherical particle system with 4 cm diameter obtained by the numerical simulation and the experimental results is shown in Figure18.Figure 18
Comparison between simulation results and experimental results of the relationship between whole load and displacement of concrete spherical particle system with 4 cm diameter.From Figure18, it can be seen that the load-displacement relationship curves of numerical simulation results and experimental results can be divided into three stages according to the increasing rate of load. In the first stage, the whole load increases slowly with the increase of displacement. In the second stage, the increasing rate of load increases, and the slope of the curve in this stage is about 1. In the third stage, the whole load increases with the increase of displacement. It increases rapidly, and the increasing rate is obviously higher than that of the first and second stages.In the first and second stages, no matter the numerical simulation results or the experimental results, with the increase of displacement, the whole load fluctuates obviously. This is due to the fact that the particles contacted with the pressure plate, they are broken up, and the fragments fill into the voids between the particles rapidly, and the whole displacement of the particle system decreases rapidly, while the falling speed of the pressure plate is less than that of the displacement of the particle system. The number of fragments in contact with the pressure plate decreases; as a result, the force on the pressure plate decreases and a short unloading process occurs. When the pressure plate contacts the particles closely again, the load begins to increase gradually. In these two stages, the error between the numerical simulation results and the experimental results is about 15%. The reason is that much powder debris is formed when the particles are broken up in the experiment. While the particles are meshed by the numerical simulation, the minimum tetrahedral element size cannot meet the requirement of the minimum debris size in the experiment. This leads to the inconsistency of the size of debris in simulation and experiment after particles break up. And the flow of debris and the filling of interstitial voids cannot be consistent. Therefore, there is a big error between the calculated whole load and the experimentally monitored whole load. In the third stage, the voids between particles have been filled by fragments, and most of the fragments have become the smallest element which cannot be broken again. At this time, the particle system becomes a material similar to the continuous medium. When the displacement increases, the whole load increases rapidly and there is no obvious fluctuation phenomenon. In this stage, the simulation results are very close to the experimental results, and the error is less than 5%.


#### 4.3.3. Fragmentation Size Distribution

In the numerical simulation, the fragment size statistics algorithm is developed. By obtaining the position and size information of all tetrahedrons that make up the fragments, the shape and size of the fragments are determined, and then whether the fragments can pass through the sieve hole is judged. By calculation, after the concrete spherical particle system with 4 cm diameter is broken, the smallest fragment passes through a 6 mm aperture sieve. By calculating the volume of the fragments of each sieve in turn and dividing by the sum of the volume of all the particles in the particle system, the percentage of the mass of the fragments smaller than a certain particle size in the total mass can be obtained; thus, the size distribution of the fragments can be obtained.In the experiment, after loading, the particle system was divided into three layers according to the height. Then, take out layer by layer and conduct sieving experiment. The sieve pore size was 1 mm, 1.7 mm, 2.5 mm, 3.1 mm, 5 mm, 7.1 mm, 10 mm, 16 mm, 20 mm, 25 mm, 35.5 mm, 40 mm, and 45 mm, respectively. In the experiment, lots of powder debris were formed after the particles were crushed; they can pass through a sieve with a diameter of 1 mm, but the size of tetrahedral element in numerical simulation cannot be too small; otherwise, too many elements will affect the computational efficiency. And the size of the smallest debris in numerical simulation cannot reach the size of the smallest debris in the experiment.The comparison between the layered sieving curve of fragments after the loading of concrete spherical particle system with a diameter of 4 cm obtained by numerical simulation and the experimental results is shown in Figure19.Figure 19
Comparison of simulation results and experimental results of stratified sieving curve of concrete spherical particle system with diameter of 4 cm after loading.From Figure19, it can be seen that the size of the largest fragments produced by upper and lower particles is smaller than that of the largest fragments in the middle. The percentage of the mass of the smaller fragments in the upper and lower layers is obviously higher than that of the middle fragments. The numerical simulation results and experimental results show that the fragmentation of the granular system has the characteristics of delamination failure. Due to the limitation of element size in numerical simulation, it is difficult to obtain fragments with identical sizes and shapes, and the size of the smallest fragments is also significantly larger than that of the smallest fragments in experiment, resulting in the error between simulation results and experimental results of about 23%.


## 4.1. Numerical Simulation for Flow and Fracture of Granular Materials

### 4.1.1. Numerical Sample Generation and Mesh Generation

The method for generating numerical samples of granular materials is as follows: the bottom particles are generated according to particle size and location information, and the particles are added one by one in a certain order. Each time the particles are added, the particles will fall naturally under the acceleration of gravity. Then, the contact force between particles and the contact force between particles and boundary are calculated according to the contact information. Particles move under the action of unbalanced force. The acceleration, velocity, and displacement of particles are calculated according to the resultant force of particles. The positions of particles are adjusted until the particle system reaches equilibrium. The generated numerical sample is shown in Figures11(d)–11(f), and the particle meshing is shown in Figure 12.Figure 11
Comparisons of particle positions recorded in experiment of concrete spherical particle system with 4 cm diameter with those generated in numerical samples: (a–c) particle positions recorded in experiment; (d–f) particle positions in numerical samples.
(a)
(b)
(c)
(d)
(e)
(f)Figure 12
Particle mesh generation: (a) single particle; (b) multiple particles.
(a)
(b)


### 4.1.2. Process for Flow and Fracture of the Particle System

The concrete spherical particle system with diameter of 4 cm is loaded to 50 kN. The process of particle flow and fracture is shown in Figure13, in which the particle system has three layers from up to down like the situation in Figure 10, and there are 4 particles in each layer.Figure 13
Process for flow and fracture of concrete spherical particle system with 4 cm diameter.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)


### 4.1.3. Beam Failure Process

The disappearance process of beams between tetrahedral elements is shown in Figure14 when the concrete spherical particle system with diameter of 4 cm is loaded to 50 kN.Figure 14
Disappearance process of beams between tetrahedral elements of concrete spherical particle system with 4 cm diameter.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)From Figure14, it can be seen that when the pressure plate contacts the top particles, the top particles break up and the beams between tetrahedral elements disappear. As the load increases, the number of missing beams increases. When the loading is completed, most of the beams disappear and many separate tetrahedral elements are formed, which indicates that a large amount of powder debris is produced after the particles are broken. The number of beams disappeared in the top and bottom particles were obviously more than that in the middle layer. Most of the debris in the top and bottom layer contained only two or three beams, and the beam network composed of multiple beams could be clearly seen in the middle layer, which indicated that the degree of broken particles in the top and bottom layer was more serious than that in the middle layer.


## 4.1.1. Numerical Sample Generation and Mesh Generation

The method for generating numerical samples of granular materials is as follows: the bottom particles are generated according to particle size and location information, and the particles are added one by one in a certain order. Each time the particles are added, the particles will fall naturally under the acceleration of gravity. Then, the contact force between particles and the contact force between particles and boundary are calculated according to the contact information. Particles move under the action of unbalanced force. The acceleration, velocity, and displacement of particles are calculated according to the resultant force of particles. The positions of particles are adjusted until the particle system reaches equilibrium. The generated numerical sample is shown in Figures11(d)–11(f), and the particle meshing is shown in Figure 12.Figure 11
Comparisons of particle positions recorded in experiment of concrete spherical particle system with 4 cm diameter with those generated in numerical samples: (a–c) particle positions recorded in experiment; (d–f) particle positions in numerical samples.
(a)
(b)
(c)
(d)
(e)
(f)Figure 12
Particle mesh generation: (a) single particle; (b) multiple particles.
(a)
(b)


## 4.1.2. Process for Flow and Fracture of the Particle System

The concrete spherical particle system with diameter of 4 cm is loaded to 50 kN. The process of particle flow and fracture is shown in Figure13, in which the particle system has three layers from up to down like the situation in Figure 10, and there are 4 particles in each layer.Figure 13
Process for flow and fracture of concrete spherical particle system with 4 cm diameter.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)


## 4.1.3. Beam Failure Process

The disappearance process of beams between tetrahedral elements is shown in Figure14 when the concrete spherical particle system with diameter of 4 cm is loaded to 50 kN.Figure 14
Disappearance process of beams between tetrahedral elements of concrete spherical particle system with 4 cm diameter.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)From Figure14, it can be seen that when the pressure plate contacts the top particles, the top particles break up and the beams between tetrahedral elements disappear. As the load increases, the number of missing beams increases. When the loading is completed, most of the beams disappear and many separate tetrahedral elements are formed, which indicates that a large amount of powder debris is produced after the particles are broken. The number of beams disappeared in the top and bottom particles were obviously more than that in the middle layer. Most of the debris in the top and bottom layer contained only two or three beams, and the beam network composed of multiple beams could be clearly seen in the middle layer, which indicated that the degree of broken particles in the top and bottom layer was more serious than that in the middle layer.


## 4.2. Experimental Verification

In order to verify the validity of the simulation program more intuitively, experiments on flow and fracture of concrete spherical particle system are carried out under the conditions the same to the simulation.The layering position of the experimental particles and the layering and crushing of the particles after loading are shown in Figure15, and the layering sieving curve and the whole sieving curve of the particle system after loading are shown in Figure 16. The whole load-displacement curve and the local load-displacement curve of the granular system during loading are shown in Figure 17.Figure 15
Particle placement before loading and the particle stratification and fragmentation after loading in the experiment: (a–c) particle placement; (d–f) surface fragmentation of each layer; (g–i) fragmentation of each layer after taking out the layers.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)Figure 16
Layered sieving curve and overall sieving curve of particle system after experimental loading.Figure 17
Whole load-displacement relationship curve and local load-displacement relationship curve of experimental granular system: (a) whole load-displacement curve; (b) local load-displacement curve.
(a)
(b)From Figures15 and 16, it can be seen that the surface particles of 4 cm concrete spherical particle system are broken most seriously under 50 kN load. Most of the fragments are in powder shape. The mass of fragments with a size less than 1 mm accounts for 28% of the total mass of the layer, the fragments with a size less than 10 mm account for 59%, and the maximum particle size of the fragments is less than 25 mm. The fragmentation degree of the lower particles is relatively light, and larger fragments can also be seen. It is found that the mass of fragments with a size less than 1 mm accounts for 24% of the total mass of the layer, the fragments with a size less than 10 mm accounts for 41%, and the maximum particle size of the fragments exceeds 45 mm. The fragmentation degree of the intermediate layer is the lightest, and the mass of fragments with a size less than 10 mm accounts for only 5% of the total mass of the layer. The maximum particle size of the fragments exceeds 45 mm. Most of the particles remain in their original shape, some of them only produce cracks or have smaller fragments peeling off the particles, and some of them break up. The average size of the fragments is much larger than that of the upper and lower layers.In the experiment, the particle system is loaded according to the displacement of the pressure plate. The loading rate is 70 mm/min, and the loading value of the granular system can be read by pushing machine during the loading process. When the load value reaches 50 kN, the loading ends. From Figure17(a), it can be seen that the whole load-displacement relationship curve presents a zigzag shape, to rise gradually in a fluctuating manner, and the curve fluctuates greatly in the initial stage. When the displacement of the pressure plate reaches 6 cm, the fluctuation range of the curve decreases obviously, and the increasing rate of load also increases significantly. When the load value is 50 kN, the displacement of the particle system is 6.8 cm. As can be seen from Figure 17(b), the local load-displacement curves obtained from the two monitoring points are quite different in both trend and magnitude; the reason for this difference is mainly related to the gradation and arrangement of the particle system. The particle system is arranged in a random manner. Under the action of external loads, the force chains with different strength will be formed inside the particle system so that the forces on the bottom of the particle system are not uniform, and the loads are not distributed equally according to the number of contact points, which is significantly different from the forces on continuous media.


## 4.3. Comparison between the Simulation Results and the Experimental Results

The comparison between the experimental results and the simulation results includes the following three aspects: the position of particle system formation, the relationship between whole load and displacement, and the scale distribution of fragments.


### 4.3.1. Particle System Formation

The experimental photographs of particle positions in a 4 cm diameter concrete spherical particle system are compared with the generated numerical samples as shown in Figure11.From Figure11, it can be seen that the position of particles in the numerical samples of the particle system generated by the preprocessing program is basically the same as that recorded in the experiment. The particles generated by the program are standard spheres, while the concrete spheres in the experiment are manufactured by the die. The regularity of the spheres is different from that of the particles generated by the numerical simulation, but the error is not large. The generated particle system can meet the requirements of simulating granular materials’ flow and fracture under the same conditions.


### 4.3.2. Relationship between Whole Load and Displacement

During the experiment, the particle system was loaded by moving the position of the pressure plate, and the descending speed of the pressure plate was 70 mm/min, which was consistent with the numerical simulation loading rate. In the process of loading the particle system, the contact forces between the tetrahedral elements and the pressure plate are calculated separately. The total contact forces are summed up as the resultant forces of the pressure plate. Finally, the component of the resultant forces perpendicular to the direction of the pressure plate is calculated, and the whole load of the particle system can be obtained. Based on the displacement of the pressure plate and the whole load of the particle system, the relationship curve between the whole load and the displacement of the particle is drawn and compared with the experimental results. The comparison of the relationship curve between the whole load and the displacement of concrete spherical particle system with 4 cm diameter obtained by the numerical simulation and the experimental results is shown in Figure18.Figure 18
Comparison between simulation results and experimental results of the relationship between whole load and displacement of concrete spherical particle system with 4 cm diameter.From Figure18, it can be seen that the load-displacement relationship curves of numerical simulation results and experimental results can be divided into three stages according to the increasing rate of load. In the first stage, the whole load increases slowly with the increase of displacement. In the second stage, the increasing rate of load increases, and the slope of the curve in this stage is about 1. In the third stage, the whole load increases with the increase of displacement. It increases rapidly, and the increasing rate is obviously higher than that of the first and second stages.In the first and second stages, no matter the numerical simulation results or the experimental results, with the increase of displacement, the whole load fluctuates obviously. This is due to the fact that the particles contacted with the pressure plate, they are broken up, and the fragments fill into the voids between the particles rapidly, and the whole displacement of the particle system decreases rapidly, while the falling speed of the pressure plate is less than that of the displacement of the particle system. The number of fragments in contact with the pressure plate decreases; as a result, the force on the pressure plate decreases and a short unloading process occurs. When the pressure plate contacts the particles closely again, the load begins to increase gradually. In these two stages, the error between the numerical simulation results and the experimental results is about 15%. The reason is that much powder debris is formed when the particles are broken up in the experiment. While the particles are meshed by the numerical simulation, the minimum tetrahedral element size cannot meet the requirement of the minimum debris size in the experiment. This leads to the inconsistency of the size of debris in simulation and experiment after particles break up. And the flow of debris and the filling of interstitial voids cannot be consistent. Therefore, there is a big error between the calculated whole load and the experimentally monitored whole load. In the third stage, the voids between particles have been filled by fragments, and most of the fragments have become the smallest element which cannot be broken again. At this time, the particle system becomes a material similar to the continuous medium. When the displacement increases, the whole load increases rapidly and there is no obvious fluctuation phenomenon. In this stage, the simulation results are very close to the experimental results, and the error is less than 5%.


### 4.3.3. Fragmentation Size Distribution

In the numerical simulation, the fragment size statistics algorithm is developed. By obtaining the position and size information of all tetrahedrons that make up the fragments, the shape and size of the fragments are determined, and then whether the fragments can pass through the sieve hole is judged. By calculation, after the concrete spherical particle system with 4 cm diameter is broken, the smallest fragment passes through a 6 mm aperture sieve. By calculating the volume of the fragments of each sieve in turn and dividing by the sum of the volume of all the particles in the particle system, the percentage of the mass of the fragments smaller than a certain particle size in the total mass can be obtained; thus, the size distribution of the fragments can be obtained.In the experiment, after loading, the particle system was divided into three layers according to the height. Then, take out layer by layer and conduct sieving experiment. The sieve pore size was 1 mm, 1.7 mm, 2.5 mm, 3.1 mm, 5 mm, 7.1 mm, 10 mm, 16 mm, 20 mm, 25 mm, 35.5 mm, 40 mm, and 45 mm, respectively. In the experiment, lots of powder debris were formed after the particles were crushed; they can pass through a sieve with a diameter of 1 mm, but the size of tetrahedral element in numerical simulation cannot be too small; otherwise, too many elements will affect the computational efficiency. And the size of the smallest debris in numerical simulation cannot reach the size of the smallest debris in the experiment.The comparison between the layered sieving curve of fragments after the loading of concrete spherical particle system with a diameter of 4 cm obtained by numerical simulation and the experimental results is shown in Figure19.Figure 19
Comparison of simulation results and experimental results of stratified sieving curve of concrete spherical particle system with diameter of 4 cm after loading.From Figure19, it can be seen that the size of the largest fragments produced by upper and lower particles is smaller than that of the largest fragments in the middle. The percentage of the mass of the smaller fragments in the upper and lower layers is obviously higher than that of the middle fragments. The numerical simulation results and experimental results show that the fragmentation of the granular system has the characteristics of delamination failure. Due to the limitation of element size in numerical simulation, it is difficult to obtain fragments with identical sizes and shapes, and the size of the smallest fragments is also significantly larger than that of the smallest fragments in experiment, resulting in the error between simulation results and experimental results of about 23%.


## 4.3.1. Particle System Formation

The experimental photographs of particle positions in a 4 cm diameter concrete spherical particle system are compared with the generated numerical samples as shown in Figure11.From Figure11, it can be seen that the position of particles in the numerical samples of the particle system generated by the preprocessing program is basically the same as that recorded in the experiment. The particles generated by the program are standard spheres, while the concrete spheres in the experiment are manufactured by the die. The regularity of the spheres is different from that of the particles generated by the numerical simulation, but the error is not large. The generated particle system can meet the requirements of simulating granular materials’ flow and fracture under the same conditions.


## 4.3.2. Relationship between Whole Load and Displacement

During the experiment, the particle system was loaded by moving the position of the pressure plate, and the descending speed of the pressure plate was 70 mm/min, which was consistent with the numerical simulation loading rate. In the process of loading the particle system, the contact forces between the tetrahedral elements and the pressure plate are calculated separately. The total contact forces are summed up as the resultant forces of the pressure plate. Finally, the component of the resultant forces perpendicular to the direction of the pressure plate is calculated, and the whole load of the particle system can be obtained. Based on the displacement of the pressure plate and the whole load of the particle system, the relationship curve between the whole load and the displacement of the particle is drawn and compared with the experimental results. The comparison of the relationship curve between the whole load and the displacement of concrete spherical particle system with 4 cm diameter obtained by the numerical simulation and the experimental results is shown in Figure18.Figure 18
Comparison between simulation results and experimental results of the relationship between whole load and displacement of concrete spherical particle system with 4 cm diameter.From Figure18, it can be seen that the load-displacement relationship curves of numerical simulation results and experimental results can be divided into three stages according to the increasing rate of load. In the first stage, the whole load increases slowly with the increase of displacement. In the second stage, the increasing rate of load increases, and the slope of the curve in this stage is about 1. In the third stage, the whole load increases with the increase of displacement. It increases rapidly, and the increasing rate is obviously higher than that of the first and second stages.In the first and second stages, no matter the numerical simulation results or the experimental results, with the increase of displacement, the whole load fluctuates obviously. This is due to the fact that the particles contacted with the pressure plate, they are broken up, and the fragments fill into the voids between the particles rapidly, and the whole displacement of the particle system decreases rapidly, while the falling speed of the pressure plate is less than that of the displacement of the particle system. The number of fragments in contact with the pressure plate decreases; as a result, the force on the pressure plate decreases and a short unloading process occurs. When the pressure plate contacts the particles closely again, the load begins to increase gradually. In these two stages, the error between the numerical simulation results and the experimental results is about 15%. The reason is that much powder debris is formed when the particles are broken up in the experiment. While the particles are meshed by the numerical simulation, the minimum tetrahedral element size cannot meet the requirement of the minimum debris size in the experiment. This leads to the inconsistency of the size of debris in simulation and experiment after particles break up. And the flow of debris and the filling of interstitial voids cannot be consistent. Therefore, there is a big error between the calculated whole load and the experimentally monitored whole load. In the third stage, the voids between particles have been filled by fragments, and most of the fragments have become the smallest element which cannot be broken again. At this time, the particle system becomes a material similar to the continuous medium. When the displacement increases, the whole load increases rapidly and there is no obvious fluctuation phenomenon. In this stage, the simulation results are very close to the experimental results, and the error is less than 5%.


## 4.3.3. Fragmentation Size Distribution

In the numerical simulation, the fragment size statistics algorithm is developed. By obtaining the position and size information of all tetrahedrons that make up the fragments, the shape and size of the fragments are determined, and then whether the fragments can pass through the sieve hole is judged. By calculation, after the concrete spherical particle system with 4 cm diameter is broken, the smallest fragment passes through a 6 mm aperture sieve. By calculating the volume of the fragments of each sieve in turn and dividing by the sum of the volume of all the particles in the particle system, the percentage of the mass of the fragments smaller than a certain particle size in the total mass can be obtained; thus, the size distribution of the fragments can be obtained.In the experiment, after loading, the particle system was divided into three layers according to the height. Then, take out layer by layer and conduct sieving experiment. The sieve pore size was 1 mm, 1.7 mm, 2.5 mm, 3.1 mm, 5 mm, 7.1 mm, 10 mm, 16 mm, 20 mm, 25 mm, 35.5 mm, 40 mm, and 45 mm, respectively. In the experiment, lots of powder debris were formed after the particles were crushed; they can pass through a sieve with a diameter of 1 mm, but the size of tetrahedral element in numerical simulation cannot be too small; otherwise, too many elements will affect the computational efficiency. And the size of the smallest debris in numerical simulation cannot reach the size of the smallest debris in the experiment.The comparison between the layered sieving curve of fragments after the loading of concrete spherical particle system with a diameter of 4 cm obtained by numerical simulation and the experimental results is shown in Figure19.Figure 19
Comparison of simulation results and experimental results of stratified sieving curve of concrete spherical particle system with diameter of 4 cm after loading.From Figure19, it can be seen that the size of the largest fragments produced by upper and lower particles is smaller than that of the largest fragments in the middle. The percentage of the mass of the smaller fragments in the upper and lower layers is obviously higher than that of the middle fragments. The numerical simulation results and experimental results show that the fragmentation of the granular system has the characteristics of delamination failure. Due to the limitation of element size in numerical simulation, it is difficult to obtain fragments with identical sizes and shapes, and the size of the smallest fragments is also significantly larger than that of the smallest fragments in experiment, resulting in the error between simulation results and experimental results of about 23%.


## 5. Conclusions

In this paper, numerical simulation is used as the main research method to study the flow and fracture process of granular materials. Based on the existing discrete element method (DEM) theory, the particle system mesh generation algorithm, the block element contact detection algorithm, the beam element generation and action algorithm, and the fragment size statistics algorithm are developed. Combining the above algorithm with the existing numerical algorithm, a simulation program for the flow and fracture of granular materials is developed. The flow and fracture process of concrete granular system under external loading was simulated by the developed simulation program, and supplementary experiment under the numerical simulation conditions was carried out. Through comparative analysis, the numerical simulation results are in good agreement with the experimental results, which verifies the reliability of the developed numerical algorithm and simulation program.



---

*Source: 4562325-2019-08-15.xml*