Abstract
Gangue particles (GP) are an important part of solid filling materials in coal mines. The compression characteristics (CC) of gangue determine whether it can effectively control roof subsidence. The particle size distribution (PSD) of GP is the main factor affecting the CC; therefore, it is important to find the optimal size distribution of GP and to investigate the macrodeformation and micromotion characteristics of gangue compression. Here, Talbol theory was used to study the compression resistance of gangue granules. It is concluded that the compression modulus of continuously graded gangue is the largest when the Talbol coefficient n is 0.4. The engineering discrete element method was used to simulate and analyze the optimum PSD (n = 0.4) and to study the stress transfer of GP during compression. The results show that with the increase of stress, the microstructure of gangue particles changes in the support skeleton, the skeleton is destroyed and particles flow, thus forming a more stable support skeleton. The resultant force direction of particles changes from the initial vertical downward to the scattered distribution of the central axis and finally to a generally scattered distribution. The number of strong chains and weak chains increases, and the main conductive stress on strong chains becomes a uniform conductive stress on the weak chains. Most of the particles in the upper and middle parts of the model exhibit linear motion. The trajectories of the middle and lower particles in the model are clustered, undergoing only small displacement.
1. Introduction
The exploitation of coal resources is often accompanied by the production of large-scale solid waste, such as gangue and fly ash. Among them, the total mass of gangue deposited in China reached 5.5 billion tons, covering an area of 15,000 hectares [1, 2], causing serious environmental pollution and land waste; because of the limitations of the traditional collapse mining method, a large number of coal pillars [3–5] are stored under buildings, railways, and water bodies, causing huge waste of the coal resources. Gangue has the characteristics of low compressibility and high strength [6–8]. It is an excellent filling aggregate and can replace coal as a support to overlying strata and control roof subsidence. Meanwhile, as an effective mining waste treatment method, it has also achieved good results in solving some environmental problems caused by mining activities. Therefore, the study of gangue granules is important for environmental protection and improved coal resource recovery.
At present, the PSD is deemed to be one of the important factors affecting the CC of gangue. Zhang et al. [9] conducted compaction tests of coal gangue under different stresses and PSD conditions. It is concluded that the PSD of coal gangue samples with different lithology has fractal characteristics. Zhou et al. [10] used a SANS test machine and a steel cylinder to carry out compression tests on coal gangue of different gradings under different loading rates. The results show that the particle size and loading rate of coal gangue exert a significant influence on the deformation and energy dissipation of coal gangue. Guo et al. [11] studied the effect of particle composition on the creep of gangue under step loading by creep compaction testing of crushed stone with different particle compositions. Li and Ju [12] used a rock mechanics test system equipped with a self-made circular cylindrical piece of apparatus to compact granulated gangue backfilling materials with various binders (fly ash, cement, and lime). The influence of binders on the compaction characteristics of gangue material was studied. Zhang et al. [13] studied the formation and compaction characteristics of the caving zone. The results show that the reduction in block size and rearrangement of the fill are the main factors affecting the compaction process, which further affects the macromechanical properties of the caving zone. Zhao et al. [14] evaluated the controlling effect under the condition of waste rock filling mining, and in situ investigation of fracture expansion and development of overburden strata were undertaken. In addition, displacement and stress sensors used for monitoring roof subsidence and filling body stress were arranged behind the longwall face to monitor and record the force and deformation of such gangue. Guo et al. [15] studied compression tests of two kinds of gangue, finding that a logarithmic function is suitable for describing its equal-time stress-strain relationship, and modified the Singh–Mitchell creep model. Ma et al. [16] analyzed the effects of initial PSD and water content on the breaking behavior and permeability of granular gangue. Huang et al. [17] used PFC3D to study the macromechanical response and micromovement characteristics of loose gangue under different PSDs, confining pressures and different loading rates. Luo et al. [18] used a Hopkinson pressure bar to elucidate the effects of initial particle size and water content on sand volume and skewness behavior at high stress and high strain rates. Ma et al. [19] adopted the MTS815.02 system and a self-designed water-flow apparatus to demonstrate the influence of PSD on the compaction and seepage behaviors of a crushed limestone particle mixture. Bagherzadeh-Khalkhali and Mirghasemi [20] studied the influence of particle size on shear strength of coarse-grained soils using experimental tests in different scales and numerical simulations based on the discrete element method (DEM). Ma et al. [21] used the self-improved test system to conduct seepage-induced erosion tests under different initial porosities, particle size gradations, and water pressures. Despite many useful findings, most of abovementioned research programs focused on the same particle shape rather than different shapes of particles. Although many advances have been made in existing research, but the research on the macrodeformation and micromotion characteristics of the optimal size distribution of different shapes of GP remains sparse.
In this paper, the CC of gangue particles with the PSD matching the Talbol formula were studied through physical experiments on the basis of previous work, and the optimum PSD of GP was found. We investigated the macro- and microcompaction characteristics of the optimal PSD, through the combination of physical experiment, an engineering discrete element method (EDEM) and the Fourier series method; the shape of GP was classified to improve the calculation accuracy of EDEM. The GP under compression were investigated, and the transfer of stress, the mutual transformation of force chains, and the micromotion characteristics were explored. The compression and deformation characteristics of GP were analyzed from a microperspective, and the relationship between the development of the force chain and the macromechanical response was established, which can provide a reference for the study of the physicomechanical properties of GP.
2. Principle of Discrete Element Method and Theory of Granular Size Distribution
The principle of the discrete element method [22, 23] is to divide a complex dense particle set into independent units and separately study the interaction between each particle unit using Newton’s laws of kinematics and update each unit using a looped iterative method. The position of the particles allows us to observe the motion state of the entire particle set. According to the relationship between force and displacement, the motion displacement of the particles and the force state of the particles were obtained by Newtonian kinematics:where is the acceleration of particle i (m/s2); is the angular acceleration of particle i (rad/s2); is the mass of particle i (kg); is the moment of inertia of particle i (kg·m2); is the resultant force on particle i (N); and is the moment about the centroid of particle i (kg·m2).
The central difference method is used to integrate formula (1), and the update speed is expressed by the midpoint of the iteration time step, so as to calculate the displacement in the next time step:where N is the corresponding time(s) and is the time step.
The contact force required can be obtained by introducing the results of formula (2) into the contact model. After several iterations, the movement state of particles in a certain period of time can be analyzed.
Granular media are in a state between liquid and solid, and the movement of particles is certain, but when there are more small particles, coarse particles can only be dispersed between the smaller particles and cannot form a supporting structure; when there are more coarse particles and fewer small particles, small particles cannot fill the pore of large particles, and the skeleton formed by coarse particles is unstable. Therefore, the Talbol formula is used to guide the distribution of gangue. Fuller’s theory [24] and Schlangen [25] proposed that the closer the PSD curve is to a parabola, the greater the density and the smaller the void ratio. Talbol modified Fuller’s theory and proposed the Talbol formula:where P is percentage passing a certain sieve aperture size of different particles of gangue; d is the particle size of gangue at different levels (mm); D is the maximum particle size of gangue (mm); and n is an experimental coefficient (0.3 <n <0.7, as shown in Table 1), which is a parameter to characterize the PSD of gangue. The increase of n means that the GP with a small particle size are relatively reduced. The size of n is related to the porosity, density, and arrangement of particles, but not to the physical and chemical characteristics of the particles.
3. Experimental Work
3.1. Test Method
Using sieving distributions of 0–10 mm, 10–20 mm, 20–30 mm, and 30–40 mm (via vibrating screens), the gangue from a coal mine in Shanxi Province was screened and graded. In the heading face, coring for intact rock allowed production of standard specimens and grinding the surface of the gangue specimen ensured accuracy therein. The ratio of GP was calculated by using the Talbol formula. The grading of the gangue with n values of 0.3, 0.4, 0.5, 0.6, and 0.7 (Table 1) was used to place 2 kg gangue into the steel drum (diameter, 120 mm; height, 200 mm) for compression testing. The standard gangue specimens and the loading test of gangue granules under distribution by Talbol formula was carried out by using a Changchun Xinke YA-600 automatic pressure tester controlled by a computer. The displacement control mode was used in the test. The loading rate was 0.001 m/s. The normal load was provided by a hydraulic jack. The mechanical properties of the standard gangue specimens and the compression properties of the gangue granular media under complete lateral restraint were measured, respectively.
3.2. Compression Tests on Gangue under Complete Lateral Restraint
Figure 1 shows the stress-strain curves of gangue CC under Talbol formula distribution in laboratory tests. All graded gangue aggregates flow to form a supporting skeleton that is then destroyed with particle breakage and then flow, leading to a more stable supporting skeleton:(1)When the stress is 0-1 MPa, the strain and stress show a linear relationship, and the strain changes rapidly. At this stage, the pores between the gangues are larger and there is no bearing structure between the particles. Under low external load, the small gangues flow to the pore, and the large gangues rearrange and gradually occlude each other to form a support skeleton.(2)When the stress is 1–4 MPa, the strain rate of the bulk medium decreases as the stress increases, and the particles occlude each other to form a relatively stable supporting skeleton; when the bearing capacity of the supporting skeleton formed by the biting of the GP reaches its limit, the GP dislocate each other or the edges and corners of the particles in the gangue biting structure are broken, and the old skeleton is destroyed. The flow combination of large granular gangue forms a new skeleton structure. The broken small granular gangue supports the gangue skeleton, and the new supporting skeleton ends up with a greater bearing capacity.(3)When the stress is 4–12 MPa, the strain increases very slowly. The gangue granules have formed a stable skeleton, and the pores are filled with small particles; however, some of the gangue granules reach their ultimate strength and undergo deformation and destruction.(4)The effect of strain on the PSD is different under the same stress because of the different distribution of particles. The main reason is that the gap between the large particles is filled by the small particles, with few pores, close relationship between particles, and high stiffness. Among them, n = 0.4 distribution is reasonable, its strain is the smallest, and its stiffness is the greatest. The distribution at n = 0.3 contains more small particles; the number of large particles is relatively small and the connection between the large particles is weak and cannot form an effective support structure. The small particles are less numerous in the distribution at n = 0.6 and cannot fill the pores between large particles, so the stiffness is low.(5)There are more large particles in the distribution at n = 0.7, which forms a certain supporting skeleton before loading, and the pores between skeletons are larger. With increasing load, the large particles of gangue rotate, flow, and crush sufficiently, and the number of broken particles grows; these fill the pores, so its compressive stiffness is enhanced.

4. Numerical Simulation of Compression Testing of Gangue
4.1. Establishing a Gangue Model
Scholars considered that the shape of particles was one of the most important factors affecting the macromechanical properties of particles [26–30]. The shape of GP has a large effect on the rotation and movement of GP and increases the compaction resistance. Consequently, classification of GP is mainly based on the quantitative mathematical expression of geometric parameters such as length, width, and height. Mathematical expression is often analyzed by Fourier series. The results are displayed by the shape index and coefficient. In this experiment, the gangue was screened using a screening mesh with different particle sizes. The representative GP were measured and photographed. The shape and quantity of the GP in different sections were counted. The results show that the GP mainly have five shapes, specifically being conical, prismatic, ellipsoidal, sheet, or cylindrical in shape. Image analysis of typical GP was carried out, and the curve function of the shape of GP was obtained by fitting. ProE software was used to stack several groups of closed curves into curved surfaces, and finally a three-dimensional model of GP was established. EDEM software was used to import the three-dimensional model of GP. The spherical particles which are internally tangent to the model were placed in turn until the combined shape of all the spheres is close to the required typical particle shape. The shape diagram of typical GP and the corresponding three-dimensional model diagram are shown in Figure 2.

(a)

(b)

(c)

(d)

(e)

(f)
A Hertz–Mindlin (no-slip) model [31–32] is used in this numerical simulation experiment; both normal and tangential forces have damping components, and the damping coefficient has a certain relationship with the impact recovery coefficient. Tangential friction follows a Coulomb friction law. Rolling friction is realized by a contact independent directional constant torque model. The model calculation formula is as follows:
4.1.1. Normal Force
where is the equivalent Young’s modulus (MPa); R∗ is the equivalent radius (mm); and δn is the normal overlap (mm2).
4.1.2. Tangential Force
where St is the tangential stiffness (N/mm2); G∗ is the equivalent shear modulus (Mpa); and δt is the tangential overlap (mm2).
4.1.3. Rolling Friction
where μr is the rolling friction coefficient; Ri is the distance from the contact point to the center of mass (mm); and ωi is the unit angular velocity vector at point of contact (rad/s).
In laboratory tests, the intrinsic parameters of gangue could be obtained by the loading test of standard gangue specimens. The intrinsic parameters of steel were selected from the EDEM material database. The intrinsic parameters of gangue and steel are listed in Table 2.
It can be seen from the aforementioned that the contact parameters between particles determine the microbehavior of sliding, contact, and occlusion between particles; because the object of numerical simulation is the bulk gangue, the contact parameters between particles are difficult to obtain directly. Here, the idea of reverse modelling was adopted to determine the microparameters of gangue. The microparameters of the numerical model were debugged by inversion several times until the angle of the numerical simulation matched the experimental results. The calibration results of the main microparameters of the numerical model of gangue granules are listed in Table 3.
4.2. Numerical Simulation Results
In EDEM, a model of the aforementioned steel drum with a diameter of 120 mm and height of 200 mm was established (this is the same size as the self-made seamless steel drum), and a circular plate with a diameter of 60 mm was created for use as a simulated indenter. We created a particle factory above the steel drum to generate 2 kg of particles in the steel drum; the PSD follows the Talbol formula (for n = 0.4) and shape distribution follows the actual statistical results. The actual statistical results of the shape distribution are as follows: the ratio of the number of ellipsoidal particles; the number of conical particles; the ratio of the number of sheet particles; the number of prismatic particles; and the ratio of the number of cylindrical particles is 3 : 1.2 : 2.3 : 1.5 : 2. The numerical simulation experiment model and the physical experiment model are illustrated in Figure 3.

(a)

(b)
The relative error between the physical test results and simulation results can be expressed as follows:where ER is the relative error; Ea is the absolute error; T is the true value; εn is the physical test result; and εp is a simulated result.
When a compressive stress of 4.5 MPa was applied, we stopped the numerical simulation calculation, plotted the stress-strain curve arising therefrom, and extracted the stress-strain curve from physical experimental data (where the Talbol coefficient n is 0.4) as the stress increases from zero to 4.5 MPa for comparative analysis. Ultimately, the fitted comparisons between the physical experiments and numerical simulations are shown in Figure 4. It can be seen from the figure that the maximum error value is 0.18 and the average value is 0.11. This shows that the method of establishing a particle model and generating particles can improve the accuracy of numerical simulation.

5. Micromechanical Analysis Based on Numerical Simulation Results
Figure 4 shows a stress-strain diagram for numerical simulation of a compression test under complete lateral restraint. It can be seen that the slope of stress-strain curve at typical points A (σ = 0.1 MPa, ε = 0.046), B (σ = 1.1 MPa, ε = 0.078), C (σ = 1.78 MPa, ε = 0.122), and D (σ = 3.7 MPa, ε = 0.173) has changed. These macrochanges are caused by the change of internal microstructure, so the internal particles in bulk are changed; therefore, the formation and transformation of contact forces and strength chains were analyzed.
5.1. Analysis of Changes in the Forces on Particles
The change of the stress state of GP reflects the change of force in the process of particle compaction. The energy transfer process in the deformation process of dispersion is described. Based on the postanalysis module of EDEM, the particle model is simplified to the vector arrow of the resultant force on the particle. The force analysis graph and vector graph at typical points A, B, C, and D are shown in Figure 5.(1)In stage O-A, the stress increases from 0 MPa to 0.1 MPa and the strain increases to 0.046. From the typical Point A of the force on the particles, it can be seen that the force on the upper left and the middle part of the numerical model is the largest at 1290 N. The support skeleton has been formed between the model particles and is mainly distributed in the upper and middle parts of the model. The supporting skeleton mainly acts to transfer stress, while other particles only flow and rearrange themselves under the stress. The direction of resultant force varies.(2)In stage A-B, the stress increases from 0.1 MPa to 1.1 MPa, and the strain increases slowly from 0.046 to 0.078. From the typical Point B of the force on the particles, it can be seen that the supporting skeleton of the upper part of the numerical model has vanished, and the maximum force on the middle part of the model is 4790 N. In the middle of the model, the supporting skeleton is formed, and the resultant force on most particles is downward, and the small particles also begin to transfer stress.(3)In stage B-C, the stress increases from 1.1 MPa to 1.78 MPa, and the strain increases rapidly from 0.078 to 0.122. From typical Point C of the force distribution on the particles, it can be concluded that the supporting skeleton of the upper particles of the model disappears completely and the load is uniformly transmitted. More particles in the middle of the model began to form a supporting framework. The overall supporting framework of the model has been formed, and the maximum force is 8610 N. The resultant force of particles around the model is downward, and the resultant force on particles inside the model changes from vertical to lateral.(4)In stage C-D, the stress increases from 1.78 MPa to 3.7 MPa, and the strain increases from 0.122 to 0.173. From the typical Point D, it can be seen that the upper and lower parts of the model bear less stress, and a dense support skeleton is formed in the middle part, with the maximum stress increasing sharply to 23,600 N. The force direction of the particles forming the supporting skeleton is extremely scattered, and the resultant force of the particles with smaller force is mostly directed to the outside, and the particles with a downwards resultant force direction are rare, so the load-carrying capacity of the whole model reaches its limit state.

(a)

(b)

(c)

(d)
5.2. Distribution of Strength-Weakness Force Chains and Their Intertransformation
The force chain is a bridge to link the macroscopic and microscopic characteristics of granular particles. Since the contact force transfer inside granular particles is nonuniform, the evolution of the force chain is of great significance to the study of GP. To explore the change of contact force between GP under compression, the contact force chain network of GP under typical stress conditions A, B, C, and D was derived in EDEM. The transformation between strong and weak force chains was analyzed and the force chain distribution at each typical point was compared. Extracting the maximum contact force of the corresponding stage, a weak force chain is defined when the contact force is less than one-third of the maximum contact force. A strong chain is defined when the contact force is greater than one-third of the maximum contact force. The red chains are connected using the solid curves to vividly describe the process of connecting the strong chains to form the overall supporting structure. The contact force chain network between typical GP is illustrated in Figure 6.

(a)

(b)

(c)

(d)
From Figure 6, it can be seen that with the increase of stress, the number of strength chains between GP increases, and the strong force chains gradually shift from the upper and lower ends of the model to the middle and lower parts, forming a stable supporting structure. The strong chain in the upper part of the model plays a role in transferring larger loads when the stress is small, but with the gradual increase in stress, the upper part of the granular particles flow and reorganize, and a stable supporting structure cannot be formed between the granules. The strong chain is transformed into a weak chain. It can be seen from the comparative analysis of Figures 4 and 6 that, with a gradual increase in stress, the strong chains are gradually connected to form a stable supporting structure. At this time, the compressive capacity of materials is gradually enhanced, and the stress-strain curves gradually flatten. It can be seen from the aforementioned phenomena that the number of the strong force chains increases gradually with the stress and the number of the strong chains increases therewith to develop from a scattered distribution to mutual interconnection, forming an integral supporting structure. At this time, the gangue is compacted, and then the axial strain does not change with the increase of external applied load. Thus, the relationship between the evolution of the microscale force chain structure and the macromechanical response (gradual compaction of materials) is established.
The number of weak force chains increases with increasing stress, and they are mainly distributed in the upper part of the granular model at the beginning. The distribution of weak chains around the strong chains in the lower part of the model is sparse, but the number of weak chains around the strong chains increases gradually with increasing stress. When the stress is small, there are many voids around the support skeleton represented by the strong chain, and the small particles have no effective contact with the support skeleton in the voids; however, when the stress increases, the voids between the support skeletons decrease in size, and the small particles form effective contact with the support skeleton, forming a weak chain, thus increasing the number of weak chains in the middle and lower parts. The total number of contact force chains and the maximum contact force inside the model are shown in Table 4.
5.3. Displacement Trajectory and Angular Velocity
The angular velocity and movement trajectory indicate the movement state of GP, which has important research significance for the deformation and instability of loose particles. As seen in Figure 7, the strain in stage O-A is mainly caused by the downward movement of particles in the upper part of the model. In the middle part of the model, the particles also move downward, albeit only slightly.

(a)

(b)

(c)

(d)

(e)

(f)
In stage A-C, the particles in the upper part of the model move downward. The middle particles move up or down. In the lower part of the model, the particles move irregularly, and their trajectories are clustered.
In stage C-D, in the sparse part of the particle movement, the particles mainly present downward linear motion, while in the dense part of the particle movement, the particle trajectories are clustered.
In the displacement trajectory, the acceleration of particles is marked. The marking shows that all particles have a certain angular velocity. The maximum angular velocity occurs when the upper particles move downward, and the movement of particles is active in this stage. The strain is mainly caused by the downward displacement of the upper particles and the gradual compaction of the middle particles. Most of the particles in the upper and middle parts of the model exhibit linear motion (see magnified view of Location A in Figure 7(e)). The trajectories of the middle and lower particles in the model are clustered (see magnified view of Location B in Figure 7(f)), with small displacement, small angular velocity, and little active motion. The GP in the upper part of the model move in a linear fashion with the largest motion range of 37.8 mm, while the particles in the lower part of the model move in a cluster-like fashion with the smallest motion range of 4 mm.
6. Conclusion
(1)Fourier series analysis was used to classify the shape of gangue, including five typical shapes: cone, flake, prism, ellipsoid, and cylinder. The number of particles with different shapes was counted. When the EDEM particle factory produces particles, their shape distribution and PSD follow the actual statistical results. The method of establishing the particle model and generating particles effectively improves the accuracy of numerical simulation.(2)The stress-strain curves of five particle sizes of GP show a three-stage pattern during compression, corresponding to the variation of internal particles as follows: particle flow arrangement-forming support skeleton-skeleton destruction and particle breakage-flow arrangement-new support skeleton.(3)Among the five PSD, the interaction between particles with Talbol formula n = 0.4 is the best, the compression performance is the best, and the compression modulus is the largest. This can be used as a guide for filling mining.(4)With the increase of stress in numerical simulation, stress concentrations gradually shift from the upper and middle parts of the model to the middle and lower parts. The resultant force direction of the particles in the middle axis of the model gradually changes from the axial direction to the radial direction, while the resultant force direction of the peripheral particles gradually changes from vertically downwards to a scattered distribution.(5)In the upper part of the numerical model, the strong chain gradually transforms into the weak chain with increasing stress, from the main conductive stress of the strong chain to the uniform conductive stress of the weak chain. In the numerical model, the number of strong chains in the lower part increases gradually and their distribution is more focused. Moreover, the pore around the strong chain decreases, and the small particles in the pore play a supporting role in the skeleton, forming a large number of weak chains.(6)The strong force chain gradually develops from a scattered distribution to mutual interconnection, forming an integral supporting structure. At this time, the gangue is compacted, and then the axial strain does not change with the increase of external applied load. Thus, the relationship between the evolution of the force chain structure and the macromechanical response was established.(7)Under load, the angular velocity of particles in the upper part of the numerical model is larger, the trajectory of particles is quasilinear, and the movement of particles is more active; the angular velocity of particles in the middle and lower parts of the model is smaller, the trajectory of particles is clustered, and the movement of particles is diminished.
Data Availability
All data and models generated or used during the study are available within the article. Some data and models referred to the previous research by Fuller and Thompson [24] and Schlangen and van Mier [25] which have been cited in the text.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This research was supported by the Key Research Development Project of Hebei Province (Grant no. 18273815D) and funded by the Research Fund of State and Local Joint Engineering Laboratory for Gas Drainage and Ground Control of Deep Mines (Henan Polytechnic University) (Grant no. SJF201804) and Subsidized Project of Innovative Ability Training for Graduate Students of Hebei Education Department (Grant no. CXZZSS2020083) and the Science and Technology Research Project of Hebei Province (Grant no. QN2017031). The authors also would like to thank the editor and the anonymous reviewers for their valuable comments, which have greatly improved this paper.