Numerical Simulation Study on Flow and Fracture of Granular Materials

Liu, Mingqing; Liu, Jun; Liu, Guangkun; Shan, Peng; Xiao, Zhimin; Zhao, Futian

doi:https://doi.org/10.1155/2019/4562325

Mathematical Problems in Engineering

On this page

Abstract Introduction Numerical Simulation Conclusions Data Availability Conflicts of Interest Acknowledgments References Copyright Related Articles

Research Article | Open Access

Volume 2019 | Article ID 4562325 | https://doi.org/10.1155/2019/4562325

Numerical Simulation Study on Flow and Fracture of Granular Materials

Mingqing Liu,¹Jun Liu,²Guangkun Liu,¹Peng Shan,¹Zhimin Xiao,¹and Futian Zhao¹

Academic Editor: Gaetano Giunta

Received06 Jun 2019

Revised26 Jul 2019

Accepted30 Jul 2019

Published15 Aug 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.

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.

(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.

(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 tetrahedron i on the three coordinate axes from the storage catalogue.
	(b) Cycle all remaining tetrahedrons after removing the tetrahedron i (for (j = 0; j < num-1; i++)).
	① Generate a directory.
	② Judge whether the projection on three coordinate axes of the boundary box of tetrahedron j 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.

(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 vertices and . of A and B are obtained.
(4)	Calculate the direction vectors of the connecting line of and ; the direction vector is stored in vecN as the unit normal vector of the common surface.
(5)	The faces sharing vertices and 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 sharing vertex and 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 vertex , calculate its direction vector, and then calculate the direction vector of the line between and ; calculate their vector product, then take it to vector product with the direction vector of the edge passing through 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 , and then repeat the above operations for all edges passing through vertex .
(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 as and ; if and , 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 Algorithms 4 and 5, respectively.

(1)	Obtain the detailed contact information of tetrahedron according to the algorithm in Section 2.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 tetrahedron j (for (j = 0; j < neighborNum; j++) (neighborNum is the number of tetrahedrons in contact with tetrahedron i)).
	If there is no beam between tetrahedron i and tetrahedron j.
	① Get the centroid of tetrahedron i 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.

(1)	The list of beams between all block elements is obtained from the algorithm in Section 2.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 beam i.
	(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 Algorithm 6.

(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 tetrahedron i.
	① If included.
	Continue with the next tetrahedron cycle.
	② Otherwise, judge whether there is a beam between tetrahedron i and other tetrahedrons.
	(i) If existed.
	(1) Record the tetrahedron number and continue to determine whether the tetrahedron is associated with other tetrahedrons (excluding tetrahedron i), 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 sieve j.
	(c) Calculate the percentage of the mass of fragments smaller than the sieve j 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 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 load-displacement curve and the local load-displacement 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 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 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 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 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. KYCX18-0567).

References

X. Bian, Y. J. Cui, and X. Z. Li, “Voids Effect on the swelling behaviour of compacted bentonite,” Géotechnique, vol. 69, no. 7, pp. 593–605, 2019.
View at: Publisher Site | Google Scholar
N. Zhang, Y. Zhang, Y. Gao, R. Y. S. Ronald, Y. X. Wu, and F. Zhang, “An exact solution for SH-wave scattering by a radially multilayered inhomogeneous semicylindrical canyon,” Geophysical Journal International, vol. 217, no. 2, pp. 1232–1260, 2019.
View at: Publisher Site | Google Scholar
P. G. D. Gennes, “Reflections on the mechanics of granular matter,” Physica A: Statistical Mechanics and its Applications, vol. 261, no. 3-4, pp. 267–293, 1998.
View at: Publisher Site | Google Scholar
P. G. D. Gennes, “Granular matter: a tentative view,” Review of Modern Physics, vol. 71, no. 2, pp. S374–S382, 1999.
View at: Publisher Site | Google Scholar
J. Tejchman, “Effect of boundary, shear rate and grain crushing on shear localization in granular materials within micro-polar hypoplasticity,” in Computer Methods in Mechanics, Springer, Heidelberg, Germany, 2010.
View at: Google Scholar
L. Mishnaevsky, K. Derrien, and D. Baptiste, “Effect of microstructure of particle reinforced composites on the damage evolution: probabilistic and numerical analysis,” Composites Science and Technology, vol. 64, no. 12, pp. 1805–1818, 2004.
View at: Publisher Site | Google Scholar
W. Hu, Z. Yin, C. Dano, and P.-Y. Hicher, “A constitutive model for granular materials considering grain breakage,” Science China Technological Sciences, vol. 54, no. 8, pp. 2188–2196, 2011.
View at: Publisher Site | Google Scholar
C. A. Tang, H. Y. Liu, W. C. Zhu et al., “Numerical approach to particle breakage under different loading conditions,” Powder Technology, vol. 143-144, pp. 130–143, 2004.
View at: Publisher Site | Google Scholar
M. Khanal, W. Schubert Jr., and J. Tomas, “Ball impact and crack propagation ? Simulations of particle compound material,” Granular Matter, vol. 5, no. 4, pp. 177–184, 2004.
View at: Publisher Site | Google Scholar
M. Khanal, W. Schubert, and J. Tomas, “Crack analysis in single plate stressing of particle compounds,” in Proceedings of the The International Conference on Computational Science—ICCS 2004, vol. 3037, pp. 364–371, Springer, Kraków, Poland, June 2004.
View at: Publisher Site | Google Scholar
M. Khanal, W. Schubert, and J. Tomas, “DEM simulation of diametrical compression test on particle compounds,” Granular Matter, vol. 7, no. 2-3, pp. 83–90, 2005.
View at: Publisher Site | Google Scholar
T. Chen, Q. Fang, Z. Wang, and W. Zhu, “Numerical simulation of compression breakage of spherical particle,” Chemical Engineering Science, vol. 173, pp. 443–454, 2017.
View at: Publisher Site | Google Scholar
M. J. Metzger and B. J. Glasser, “Numerical investigation of the breakage of bonded agglomerates during impact,” Powder Technology, vol. 217, pp. 304–314, 2012.
View at: Publisher Site | Google Scholar
C. L. Bouteiller and M. Naaim, “Aggregate breakage under dynamic loading,” Granular Matter, vol. 13, no. 4, pp. 385–393, 2011.
View at: Publisher Site | Google Scholar
K. Zheng, C. Du, J. Li, B. Qiu, L. Fu, and J. Dong, “Numerical simulation of the impact-breakage behavior of non-spherical agglomerates,” Powder Technology, vol. 286, pp. 582–591, 2015.
View at: Publisher Site | Google Scholar
A. Refahi, J. A. Mohandesi, and B. Rezai, “Discrete element modeling for predicting breakage behavior and fracture energy of a single particle in a jaw crusher,” International Journal of Mineral Processing, vol. 94, no. 1-2, pp. 83–91, 2010.
View at: Publisher Site | Google Scholar
F. K. Wittel, “Single particle fragmentation in ultrasound assisted impact comminution,” Granular Matter, vol. 12, no. 4, pp. 447–455, 2010.
View at: Publisher Site | Google Scholar
S. A. Galindo-Torres, D. M. Pedroso, D. J. Williams, and L. Li, “Breaking processes in three-dimensional bonded granular materials with general shapes,” Computer Physics Communications, vol. 183, no. 2, pp. 266–277, 2012.
View at: Publisher Site | Google Scholar
K. Zhang, S. Zhang, J. Teng et al., “3D numerical simulation of particle breakage using discrete element method,” Rock and Soil Mechanics, vol. 38, no. 7, pp. 2119–2127, 2017, in Chinese.
View at: Google Scholar
H. Peron, J. Y. Delenne, L. Laloui, and M. S. E. Youssoufi, “Discrete element modelling of drying shrinkage and cracking of soils,” Computers and Geotechnics, vol. 36, no. 1-2, pp. 61–69, 2009.
View at: Publisher Site | Google Scholar
Z. B. Tong, R. Y. Yang, A. B. Yu, S. Adi, and H. K. Chan, “Numerical modelling of the breakage of loose agglomerates of fine particles,” Powder Technology, vol. 196, no. 2, pp. 213–221, 2009.
View at: Publisher Site | Google Scholar
M. Khanal and J. Tomas, “Effect of target walls on the breakage of composite particles,” Chemical Engineering and Technology, vol. 31, no. 4, pp. 609–614, 2008.
View at: Publisher Site | Google Scholar
D. O. Potyondy and P. A. Cundall, “A bonded-particle model for rock,” International Journal of Rock Mechanics and Mining Sciences, vol. 41, no. 8, pp. 1329–1364, 2004.
View at: Publisher Site | Google Scholar
Y. Yang and S. G. Tang, “Discrete element method for particle flow simulation and its progress,” Chinese Powder Technology, vol. 5, pp. 38–43, 2006, in Chinese.
View at: Google Scholar
O. D. L. Strack and P. A. Cundall, The Distinct Element Method as a Tool for Research in Granular Media: Part I. Report to the National Science Foundation, University of Minnesota, Minneapolis, MN, USA, 1978.
P. A. Cundall and O. D. L. Strack, The Distinct Element Method as a Tool for Research in Granular Media: Part II. Report to the National Science Foundation, University of Minnesota, Minneapolis, MN, USA, 1979.
P. A. Cundall and O. D. L. Strack, “Discussion: a discrete numerical model for granular assemblies,” Géotechnique, vol. 30, no. 3, pp. 331–336, 1980.
View at: Publisher Site | Google Scholar
J. Liu, B. Yun, and C. B. Zhao, “An improved specimen generation method for DEM based on local Delaunay tessellation and distance function,” International Journal for Numerical and Analytical Methods in Geomechanics, vol. 36, no. 5, pp. 653–674, 2012.
View at: Publisher Site | Google Scholar
P. A. Cundall, “Formulation of a three-dimensional distinct element model—part I: a scheme to detect and represent contacts in a system composed of many polyhedral blocks,” International Journal of Rock Mechanics and Mining Sciences and Geomechanics Abstracts, vol. 25, no. 3, pp. 117–125, 1988.
View at: Publisher Site | Google Scholar
R. Hart, P. A. Cundall, and J. Lemos, “Formulation of a three-dimensional distinct element model-Part II: mechanical calculations for motion and interaction of a system composed of many polyhedral blocks,” International Journal of Rock Mechanics and Mining Sciences and Geomechanics Abstracts, vol. 25, no. 3, pp. 117–125, 1988.
View at: Publisher Site | Google Scholar
J. Liu, C. B. Zhao, and B. Yun, “Numerical study on explosion-induced fractures of reinforced concrete structure by beam-particle model,” Science China Technological Sciences, vol. 54, no. 2, pp. 412–419, 2011.
View at: Publisher Site | Google Scholar
J. Liu, B. Yun, and C. B. Zhao, “A grid generation algorithm for failure simulation of brittle materials,” Journal of Computational Mechanics, vol. 29, no. 2, pp. 230–235, 2012, in Chinese.
View at: Google Scholar
Y. S Jia, Continuous Simulation of Blasting Collapse Process in Demolition of Large Structures, Wuhan University of Technology, Wuhan Shi, China, 2010, in Chinese.
S. B. Xiao, “Quaternion method and its application,” Progress in Mechanics, vol. 23, no. 2, pp. 249–260, 1993, in Chinese.
View at: Google Scholar
H. Liu and J. Liu, “Comparison of three algorithms for spherical particle contact detection in discrete element method,” Journal of Heilongjiang University of Science and Technology, vol. 16, no. 6, 2006, in Chinese.
View at: Google Scholar
E. G. Nezami, Y. M. A. Hashash, D. Zhao, and J. Ghaboussi, “A fast contact detection algorithm for 3-D discrete element method,” Computers and Geotechnics, vol. 31, no. 7, pp. 575–587, 2004.
View at: Publisher Site | Google Scholar

Copyright

Copyright © 2019 Mingqing Liu et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

989

Downloads

852

Citations