A Laboratory and Numerical Simulation Study on Compression Characteristics of Coal Gangue Particles with Optimal Size Distribution Based on Shape Statistics

Gangue particles (GP) are an important part of solid ﬁlling materials in coal mines. The compression characteristics (CC) of gangue determine whether it can eﬀectively control roof subsidence. The particle size distribution (PSD) of GP is the main factor aﬀecting the CC; therefore, it is important to ﬁnd 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 coeﬃcient 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 ﬂow, 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 ﬁnally 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.


Introduction
e 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][4][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][7][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. erefore, 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. e 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). e 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. e 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 PFC 3D 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. e GP under compression were investigated, and the transfer of stress, the mutual transformation of force chains, and the micromotion characteristics were explored. e 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.

Principle of Discrete Element Method and
Theory of Granular Size Distribution e 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. e 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 u i → is the acceleration of particle i (m/s 2 ); θ i → is the angular acceleration of particle i (rad/s 2 ); m i is the mass of particle i (kg); I i is the moment of inertia of particle i (kg·m 2 ); F → is the resultant force on particle i (N); and M �→ is the moment about the centroid of particle i (kg·m 2 ). e 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 Δt is the time step. e 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. erefore, 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, 2 Mathematical Problems in Engineering 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. e increase of n means that the GP with a small particle size are relatively reduced. e size of n is related to the porosity, density, and arrangement of particles, but not to the physical and chemical characteristics of the particles.

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. e ratio of GP was calculated by using the Talbol formula. e 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. e 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. e displacement control mode was used in the test. e loading rate was 0.001 m/s. e normal load was provided by a hydraulic jack. e mechanical properties of the standard gangue specimens and the compression properties of the gangue granular media under complete lateral restraint were measured, respectively.

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. e flow combination of large granular gangue forms a new skeleton structure. e 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. e 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) e effect of strain on the PSD is different under the same stress because of the different distribution of particles. e 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. e 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. e 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) ere 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.

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][27][28][29][30]. e 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. e results are displayed by the shape index and coefficient. In this experiment, the gangue was screened using a screening mesh with different particle sizes. e representative GP were measured and photographed. e shape and quantity of the GP in different sections were counted. e 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. e 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. e shape diagram of typical GP and the corresponding three-dimensional model diagram are shown in Figure 2. 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. e model calculation formula is as follows: where E * is the equivalent Young's modulus (MPa); R * is the equivalent radius (mm); and δ n is the normal overlap (mm 2 ).

Tangential Force
where S t is the tangential stiffness (N/mm 2 ); G * is the equivalent shear modulus (Mpa); and δ t is the tangential overlap (mm 2 ).

Rolling Friction
where μ r is the rolling friction coefficient; R i 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. e intrinsic parameters of steel were selected from the EDEM material database. e 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. e microparameters of the numerical model were debugged by inversion several times until the angle of the numerical simulation matched the experimental results. e calibration results of the main microparameters of the numerical model of gangue granules are listed in Table 3.

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. e 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  e numerical simulation experiment model and the physical experiment model are illustrated in Figure 3. e relative error between the physical test results and simulation results can be expressed as follows: where E R is the relative error; E a 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. is shows that the method of establishing a particle model and generating particles can improve the accuracy of numerical simulation.

Analysis of Changes in the Forces on Particles.
e change of the stress state of GP reflects the change of force in the process of particle compaction. e 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. e force analysis graph and vector graph at typical points A, B, C, and D are shown in Figure 5.           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. e overall supporting framework of the model has been formed, and the maximum force is 8610 N. e 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. e 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.

Distribution of Strength-Weakness Force Chains and eir
Intertransformation. e 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. e 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 onethird of the maximum contact force. e red chains are connected using the solid curves to vividly describe the process of connecting the strong chains to form the overall supporting structure. e contact force chain network between typical GP is illustrated in Figure 6.
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. e 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. e 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. us, the relationship between the evolution of the microscale force chain structure and the macromechanical response (gradual compaction of materials) is established. e 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. e 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. e total number of contact force chains and the maximum contact force inside the model are shown in Table 4.

Displacement Trajectory and Angular Velocity.
e 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.
In stage A-C, the particles in the upper part of the model move downward. e 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. e marking shows that all particles have a certain angular velocity. e maximum angular velocity occurs when the upper particles move downward, and the

Conclusion
(1) Fourier series analysis was used to classify the shape of gangue, including five typical shapes: cone, flake, prism, ellipsoid, and cylinder. e 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. e method of establishing the particle model and generating particles effectively improves the accuracy of numerical simulation. (2) e 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 breakageflow 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. is 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. e 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) e 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. us, 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 ompson [24] and Schlangen and van Mier [25] which have been cited in the text.

Conflicts of Interest
e authors declare that they have no conflicts of interest.