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.
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 particlebeam model, parallel bonding model, etc.) among the elements. In twodimensional simulation, the common element shapes are circular and polygon, while in threedimensional 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 threedimensional 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 threedimensional block element beamparticle model [31], particle failure simulation mesh generation algorithm [32], and threedimensional 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.

2.2. Contact Detection Algorithm for Block Elements
2.2.1. FirstOrder Contact Detection Algorithm
Based on the theory of the boundary box method [35], this paper develops a firstorder contact detection algorithm, which is shown in Algorithm 2.

2.2.2. SecondOrder Contact Detection Algorithm
In this paper, the fast commonplane 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 secondorder contact detection algorithm for block elements is developed, as shown in Algorithm 3.

2.3. Generation and Action Algorithm of Beam Elements
The algorithm of beam element generation and action adopted in this paper are shown in Algorithms 4 and 5, respectively.


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 Algorithm 6.

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 Figure 1.
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 Figure 2.
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 Figure 3.
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 Figure 4.
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 Figures 5–9.
(a)
(b)
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(e)
(f)
It can be seen from Figure 5(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 Figure 6 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 Figures 7 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 Figure 9, 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 Figure 10. 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.
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 Figures 11(d)–11(f), and the particle meshing is shown in Figure 12.
(a)
(b)
(c)
(d)
(e)
(f)
(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 Figure 13, 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.
(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 Figure 14 when the concrete spherical particle system with diameter of 4 cm is loaded to 50 kN.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
From Figure 14, 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 Figure 15, and the layering sieving curve and the whole sieving curve of the particle system after loading are shown in Figure 16. The whole loaddisplacement curve and the local loaddisplacement curve of the granular system during loading are shown in Figure 17.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(a)
(b)
From Figures 15 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 Figure 17(a), it can be seen that the whole loaddisplacement 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 loaddisplacement 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 Figure 11.
From Figure 11, 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 Figure 18.
From Figure 18, it can be seen that the loaddisplacement 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 Figure 19.
From Figure 19, 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.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This study was supported by the National Natural Science Foundation of China (grant no. 51874118) and the Fundamental Research Funds for the Central Universities (no. KYCX180567).