Numerical Study of Spherical Particle Breakage Based on Nonlinear Discrete Interface Bonded Model

To study the contact mechanical characteristics and breakage mechanism of single particle, the particle contact experiment is simulated using the discrete element method. In this paper, to solve the problem that the envelope of the shear strength of the classical model is linear, a nonlinear discrete surface bonded model is introduced by applying nonlinear tangential strength criterion. The new model is compared with the classical model to verify its validity. Then, based on the new model, the particle contact experiment is numerically simulated, and the ball-plane contact breakage mechanism is analysed and discussed from the mesoscopic perspective. The results reveal that the nonlinear discrete interface bonded model can provide shear properties similar to those of actual materials. The numerical result reflects the shear failure of the edge of the conical region and the tensile failure of the whole specimen. Combined with the analysis of the results of stress distribution, velocity field, and force chain distribution, the formation mechanism of the cone core can be obtained.


Introduction
Granular materials, such as sand, gravel, ballast, and rock ll, are common in nature and are widely used in civil engineering and geotechnical engineering. e inherent characteristics of natural granular materials, such as brittleness and natural defects, cause the particles to break under external forces, especially in a state of high stress. Particle breakage will create changes in the mechanical properties of particles, which may cause a slew of engineering issues. erefore, the study on the breakage behavior and mechanism of coarse-grained soil plays a very important role in the engineering application of coarsegrained soil. However, due to the limitations of experimental methods, research on mesoscale particle breakage based on physical experiments is greatly restricted. Numerical simulation analysis provides an e ective method for studying the particle breakage characteristics at different scales.
Nowadays, many numerical methods have been introduced for particle breakage studying. Among them, the discrete element method (DEM) can establish specimen of any shape and can simulate the formation of cracks. e DEM method is computationally e ective because it does not require the generation of a complex particle grid. After many years of development, the DEM-based particle breakage simulation has been widely used, which provides an e ective way to study the mesoscopic mechanism of particle breakage. Many scholars [1][2][3][4][5][6][7][8][9][10] have studied particle breakage based on the DEM and have obtained a large number of research results.
At present, research on particle breakage based on the DEM mainly focuses on the in uence of particle breakage on the macroscopic mechanical properties of particle groups [11][12][13][14][15][16][17][18][19]. In the existing discrete -element analysis of single particle breakage, scholars mainly focus on the simulation of macroscopic breakage phenomena [20][21][22][23][24], as well as the selection and calibration of discrete element meso-parameters [25,26] and breakage criteria [27,28] in the simulation. However, there are few studies on the mechanical mechanism of single particle breakage from the meso-point of view of the stress distribution and force chain distribution based on the DEM.
However, the two current particle breakage simulation methods based on the DEM, that is, the particle breakage simulation method based on the fragment replacement method (FRM) [29] and bonded particle model (BPM), also have some limitations. e particle breakage simulation method based on the FRM can reflect the mechanical response of particle aggregates to a certain extent, which provides an effective method for large-scale physical experiments and mesoscale research [30][31][32]. However, in numerical simulations based on the FRM, the number of replacement fragments, the size of the fragments, and the replacement logic have a substantial impact on the simulation efficiency and results [33,34]. Moreover, numerical simulations based on the FRM not only require that the particle breakage criterion is obtained in advance, but problems such as the nonconservation of mass after breakage replacement need to be solved [35][36][37][38][39], which makes it difficult to meet the needs of mesoscopic analysis of the particle breakage mechanism. At present, numerical simulations based on the FRM are mainly used to study the influence of particle breakage on the overall characteristics in large particle swarm simulations. e particle breakage simulation method based on the BPM can simulate the particle breakage behaviour to a certain extent and can effectively reproduce the complex mechanical response of breakable particle aggregates [40,41].
is method can simulate complex particle shapes [42][43][44]. However, the numerical simulation of particle breakage based on the traditional BPM often has difficulty replicating the macroscopic mechanical properties of actual materials. It can only simulate the experimental phenomena, which makes it difficult to provide the basis for the analysis of the breakage mechanism [45][46][47].
In this paper, aiming at improving the deficiency of the classical discrete-element bond model, the nonlinear discrete interface bonded model is obtained by improving the tangential bond strength criterion based on the flat-joint contact model. e ball-surface particle contact experiment of a 50 mm spherical gypsum sample is simulated. e simulation results are analysed from the force-displacement curve, breaking morphology, stress distribution, and force chain distribution, and then the meso-mechanism of particle breaking is studied.

Nonlinear Discrete Interface Bonded Model
As the classical BPM in PFC, the parallel bonded model is widely used to simulate rock and rock-like materials because it can effectively simulate cement. e macroscopic phenomenon of numerical simulation of particle breakage based on the parallel bonded model is similar to that of laboratory test. However, due to the insufficient anti rotation capability of the contact, the details of meso-phenomenon of particle breakage are partially different from the actual test.
Moreover, the shear strength envelope of numerical simulation based on the parallel bonded model is linear, which is inconsistent with the actual mechanical properties of rock. erefore, there are some limitations in numerical simulation of particle breakage based on the parallel bonded model. It is difficult to explain the meso-mechanism of particle breakage. In view of the aforementioned shortcomings, based on the flat-joint model proposed by Potyondy [48], a nonlinear discrete interface bonded model suitable for simulating brittle rock breakage is established by introducing nonlinear shear line strength criterion.
e logical structure of the model is shown in Figure 1. e model basic particles interact with each other through the contact surface which is a rigid disk, as shown in Figure 1(a). e contact surface is discretized into 16 sectors along the radial and axial directions respectively as shown in Figure 1(b). Each sector adopts the same logical structure. e logical structure of the bonded state and the unbonded state is shown in Figure 1(c). For each discrete sector, the equivalent normal and shear stresses of the contact are as follows: where F (sec) n is the normal contact force, F (sec) s is the tangential contact force, and A (sec) is the area of the contact element.
When the normal or tangential bond of the sector is destroyed, it is considered that the bond of the sector has been destroyed. e normal force acting on the contact of discrete sector is tension or pressure. e normal bonded strength of contact is defined by tensile strength σ c . At the same time, the ultimate compressive strain parameter β is introduced, and the bond failure occurs when the equivalent compressive strain is equal to β. erefore, the arithmetic expression of normal compression bond strength σ c ′ in numerical calculation is as follows: where R is the contact radius and σ c is the normal bond strength. e default value of β is 1.0, which is equivalent to normal pressure failure without considering bonding. e contact shear strength of the classical bonded model is defined using the friction coefficient (unbonded state) or Moore-Coulomb criterion (bonded state). e intensity envelopes of both are linear. e triaxial test results based on the classical bonded model are quite different from the actual results. erefore, the nonlinear strength criterion of rock is introduced in the model.
It is generally believed that the quadratic parabolic Mohr strength criterion can comprehensively reflect the nonlinear strength characteristics of rock. is criterion can be used to 2 Advances in Materials Science and Engineering describe both plastic rock and the shear failure of brittle rock. e specific forms of the Mohr strength envelope include linear, quadratic parabolic, hyperbolic forms, and so on. Among them, the strength envelope of rock with weak to subhard lithology, such as sandstone and marl, can often be fitted by a quadratic parabola. In the nonlinear discrete interface bonded model, considering that the properties of the high-strength gypsum selected in the experiment [49] are similar to those of soft rock such as fine sandstone, the quadratic parabolic Mohr strength criterion is chosen as the tangential bond strength criterion of the model. e ultimate equivalent shear stress of the discrete unit can be expressed as follows: where σ c is the normal bond strength. e normal stress on the unit is positive under tension. To avoid a negative value in the root sign on the right side of the equation, tensile failure of the bond is first detected in each operation, and then the shear strength is determined. If the bond of the unit is damaged, the strength criterion shown in equation (3) fails immediately. e element follows the friction limit criterion shown in (4), which provides a certain tangential bearing capacity under compression.
where μ is the sliding friction coefficient of the bond. If the equivalent shear stress of the unit is greater than the tangential friction limit, the shear force of the unit maintains the peak value, and tangential relative sliding can occur. e improved model may form aggregates by bonding a certain number of elementary particles and simulate the formation of cracks by breaking the bonds between elementary particles.

Validation of Improved Model
Based on the new nonlinear discrete interface bonded model, a cylinder specimen with a diameter of 50 mm and a height of 100 mm is generated by referring to the material simulation library provided by Itasca. e minimum elementary particle size of the model is 2.5 mm, and the maximum elementary particle size is 5.0 mm. e average size is 3.75 mm. Tests of fundamental mechanics, such as direct tensile tests and triaxial tensile tests, are carried out to calibrate the macroscopic mechanical properties of the specimens. is can be used to verify that the macroscopic mechanical properties of the specimens in the improved model were in accordance with the design.
In the direct tensile test shown in Figure 2(a), the particles on the upper and lower surfaces of the sample move upward and downward at a constant velocity of 0.001 mm/s, respectively, to simulate the loading boundary conditions in the actual test. Figure 2(b) shows the test curve obtained by direct tensile simulation using the nonlinear bonded model and parallel bonded model. e two test curves in the figure are basically linear, and there is no loading platform before the failure of the sample. e failure characteristic is brittle failure.
In the triaxial compression test shown in Figure 3, the radial constraint of the specimen is controlled by the annular wall. e upper and lower walls of the sample are loaded inward at a constant speed of 0.001 mm/s. As shown in Figure 3(b), the peak deviatoric stress of the specimen increases with increasing radial restraint stress in the simulated triaxial compression test. Similar to actual rock or rock-like materials, when the confining pressure of the triaxial test is sufficiently large, the specimen based on the nonlinear n c t  c s c Interface n x (1) x c Advances in Materials Science and Engineering discrete interface bonded model will no longer show the characteristics of brittle failure, and the equivalent deviator stress will show a postpeak platform after reaching the peak value. When the confining pressure level is high, the equivalent deviatoric stress can continue to rise until reaching the peak value after the linear stage of the curve.
Considering that the shear strength of a contact bond is related to its tensile strength in the nonlinear discrete interface bonded model, the results of the triaxial test and direct tensile test are obtained as stress circles, and the strength envelopes of the stress circles are made by the arc tangent fitting method, as shown in Figure 4(a). e results show that the shape of the shear strength envelope of the specimen based on the nonlinear discrete interface bonded model is similar to that of the nonlinear Mohr shear strength envelope, which can be well fitted by the quadratic parabolic strength criterion. In this case, the expression of the shear strength envelope is as follows:   Advances in Materials Science and Engineering According to the setting of the model meso-parameters, the theoretical shear strength expression of the bonded state is as follows: In the formula, σ is positive in tension. e triaxial test results verify the e ectiveness of introducing the quadratic parabolic Mohr strength criterion into the improved bonded model. Comparing the results of triaxial tests, as shown in Figure 4, it is found that the specimen based on the nonlinear bonded model can provide a relatively higher triaxial compressive strength and shear strength under the same con ning pressure than that based on the parallel bonded model.
e macroscopic shear strength envelope shows that the shear strength envelope based on the improved model is close to the quadratic parabolic Mohr shear strength criterion. It is more suitable for simulating brittle soft rock or other materials suitable for the quadratic parabolic Mohr strength criterion. e nonlinear discrete interface bonded model e ectively improves the de ciency of the parallel bonded model and provides abilities for the study of the meso-mechanism of coarsegrained soil particle breakage.

Calibration and Verification of Numerical Model Parameters
In order to make the macroscopic mechanical properties of the numerical model well tted with the test material, it is necessary to calibrate and adjust the numerical meso-parameters. According to the correlation and sensitivity analysis of parameters [49] and referring to the parameter adjustment method in the PFC o cial manual, the crossdebugging simulation is carried out by interpolation analysis and iterative approximation. Finally, the mesoscopic parameters of the numerical simulation of gypsum particle breakage suitable for the nonlinear discrete interface bonded model are obtained, as shown in Table 1.
Based on the aforementioned parameters, the numerical models in Figure 5 are established for uniaxial compression test and splitting test. e results show that the numerical simulation matches the test well, as shown in Figure 6. It is worth explaining that the peak loading force of the numerical splitting test is basically consistent with the experiment results, but the secant sti ness in the initial stage is quite di erent.
rough the analysis of the indoor test process, it is considered that the curve sti ness is small because the loading conditions are not strictly met in the initial stage of the Brazilian splitting test. With the progress of loading, the loading force will gradually be uniformly distributed on the contact surface. At this time, the loading boundary condition of the sample is consistent with the experimental hypothesis, and the force-displacement curve is basically linear. e slope of the force-displacement curve of the numerical splitting test is generally consistent with the results of the experiment.
In conclusion, the meso-parameters based on the improved discrete contact surface bonded model are reasonable and can accurately t the fundamental mechanical properties of high-strength gypsum used in indoor contact tests. e discrete contact surface bonded model can be used to simulate the particle contact breakage process. e basic mechanical properties of the simulated sample material are compared with those of high-strength gypsum, as shown in Table 2.

Mechanical Characteristics.
For accurate simulation and e ective comparative analysis, the shape and size of the numerical model are consistent with those of the highstrength gypsum samples used in previous contact experiments [49]. e model generation parameters are shown in Table 3. e initial numerical model is shown in Figure 7.
e numerical analysis of the ball-plane particle contact experiment is performed. e maximum strain rate is set to 5 × 10 − 2 s − 1 , which can be converted into a maximum loading rate of 2.5 mm/s. e force-displacement curves recorded by the numerical simulation of the ball-plane

Advances in Materials Science and Engineering
Loading plate   [49], as shown in Figure 8. e force-displacement curve obtained by the numerical simulation is generally consistent with that of the experiments. e numerical simulation curve also effectively reflects the four loading stages divided according to the curve shape: (1) initial linear stage; (2) nonlinear stage with reduced contact stiffness; (3) contact stiffness linear stage; and (4) overall breakage after the peak. In the second stage, the contact stiffness decreases to a certain extent but then gradually recovers. e shape of the curve indicates that there is a nonlinear process in the contact between the elementary particles in the sample, that is, the failure of bonding or the separation of the contact. According to the loading stage of the experiment and the variation range of the curve slope, it can be confirmed that the change is caused by the local bond failure of the contact point. e through crack is formed and the model is broken into two parts. Among them, the bonded fractures are mostly tension failures, and the shear failures occur less, which is mostly concentrated in the conical area below the contact point. Figure 8(d) shows a side view of the final broken shape of the numerical model. As can be seen from the figure, the final breakage pattern of the specimen in the numerical simulation is also similar to that in the experiment, both of which break symmetrically along the section passing through the upper and lower loading points into almost complete parts.

Stress Distribution.
According to the results of the indoor ball-plane contact experiment, the final penetrating fracture surface of the specimen will pass through the central section. erefore, in the process of numerical simulation, 10 measurement balls are introduced vertically along the geometric centre of the sample (that is, the sphere centre) to record the horizontal stress distribution at the axis of the sample, as shown in Figure 9. e measurement balls are uniformly distributed along the diameter direction of the sample. e diameter of the measurement ball is 5 mm. When the contact force is 6.5137 kN, the recorded horizontal stress values of the measurement balls are used to draw the horizontal stress distribution map at different heights on the axis according to the sample height. e results are shown in Figure 9(a). e compressive stress is positive, and the tensile stress is negative. As the sample approaches overall  Material density 2650 kg/m 3 r min Minimum particle size 1.0 × 10 − 3 m r max Maximum particle size 1.6 × 10 − 3 m r avg Average particle size 1.

Advances in Materials Science and Engineering
breakage in the current state, the maximum equivalent tensile stress is close to the bond tensile strength used for the contact mechanics parameters. e maximum equivalent tensile stress occurs at the bottom third measurement ball with a value of 3.848 MPa.
Hiramatsu and Oka [50] obtained analytical solutions of the internal stress elasticity of spherical specimens under uniaxial loading in 1966 by photo elastic mechanical experiments and analytical methods. e theoretical horizontal stress distribution on the central longitudinal axis according to the theoretical results is shown in Figure 9(b). e equivalent horizontal stress distribution on the central longitudinal axis in the numerical simulation is similar to the theoretical distribution map.
e results show that the horizontal stress in the conical core is compressive and that the stress level is high before overall breakage. To some extent, this explains the shear fracture of the bond at the edge of the conical nucleus. e horizontal stress in the middle part of the longitudinal axis of the centre of the sample is mainly tensile stress, which explains why the bond fracture on the fracture surface is mainly tensile fracture when the nal overall fracture occurs. Figure 10 is the central longitudinal section of the instantaneous velocity eld of particle breakage. e velocity eld at the moment of particle breakage shows that the broken block tends to move away from the central longitudinal axis. e elementary particles in the broken block have similar movement rates and the same direction.

Velocity Field.
is indicates that the broken block is of strong integrity, and the movement trend is consistent with the overall tensile fracture characteristics. e breakage phenomenon is generally the same as the test phenomenon (Figure 8(e)). On the surface of the sample at the edge area of the loading plate, the elementary balls with different movement directions from the whole can be observed. It can be considered that the local breakage occurred at the corresponding position, resulting in the collapse of debris.

Force Chain Distribution.
e contact force chain in the numerical model is analysed to study the transfer and distribution of the loading force in the numerical specimen. e direction of the force chain is consistent with the direction of the maximum principal stress of the external load [51]. e particles near the cross-section perpendicular to the penetrating breakage surface are selected for observation. e distribution of the compression chain and tension chain during the loading process is shown in Figures 11 and 12, respectively. e blue and red line segments are used to represent the force chain transfer path, and the thickness of the line segment is proportional to the magnitude of the contact force. Owing to the different magnitudes of contact forces under tension and compression, different linewidth scale factors are used. e distribution diagram of the contact compression chain indicates that the contact pressure is primarily concentrated near the contact point at the beginning of loading. With continuing loading, the loading platform is produced.
rough the loading platform, the contact pressure is conveyed and diffused into the model, where it is uniformly distributed on the central horizontal section of the sample. e pressure chain in the conical area beneath the loading plate spreads in all directions. e angle formed by the force chain at the conical area's edge and the vertical direction is roughly 45°. is demonstrates that the balls in the corresponding position of the conical region are in a state of multidirectional compression. e equivalent horizontal stress is compressive stress, which is consistent with the results of stress distribution analysis. e distribution diagram of the tension chain shows that the maximum tension chain in the numerical sample is mainly concentrated near the central longitudinal axis. e direction of the tension chain is basically perpendicular to that of the compression chain except for a small area near the contact point. e number of force chains in the conical area below the contact point is significantly less than that in the surrounding area. e area with the highest force chain density appears at the bottom of the conical area before the overall breakage of the numerical model, as shown in Figure 12(b). At the instant of overall breakage, a large number of tension chains on the central longitudinal axis disappear, that is, the corresponding contact undergoes tensile failure and loses the tensile capacity. e distribution of tension chains at the moment of overall breakage reflects the arched effect, and a tension chain with an approximately uniform distribution appears on the outer surface of the sample. At the instant of overall breakage of the numerical sample, the contact force of the corresponding basic particles at the interruption of the force chain on the central axis is extremely unbalanced, and the disintegration will be accelerated under the combined external force. is reflects the brittle failure characteristics.

Discussion
Previous investigations have revealed that after the sample has been broken, there is a cone core beneath the contact location, as shown in Figure 13. During the experiment, the cone core penetrates into the ball particle, which results in the final breakage of the particle [49]. According to the aforementioned analysis, it can be seen that the final broken phenomenon of the model in the numerical simulation is consistent with the laboratory experiment, and the result well reflects the important feature of the cone core. e results of the numerical simulation show that a significant degree of shear failure occurs at the edge of the cone core. Combined with the aforementioned stress distribution, velocity field, and force chain distribution, the formation mechanism of the cone core can be obtained. Figure 14 shows the distribution of force chain during loading and the Advances in Materials Science and Engineering bonded fracture and velocity field at the moment of breakage in the cone core region. As shown in the figure, under the action of the loading force, a certain area below the loading plate is in the stress state of triaxial compression (see Figure 14(a)). When the contact force reaches the critical value, the shear stress on the oblique section in the region reaches the shear strength of the material and causes shear failure. e bonded failure is shown in Figure 14(c), where the yellow discs represent the bonded failures caused by tension, and the red discs represent the bond failures caused e shear failure angle is perpendicular to the direction of the compression chain at the edge of the conical zone, and the angle between the shear failure angle and the vertical direction is approximately 45°. With the accumulation of local shear failure, the cone core is gradually formed. At the moment of particle breakage, the instantaneous movement direction of the balls in the cone core is mainly downward, while the outer particles move to both sides (see Figure 14(d)). e whole ball is broken into two halves.

Conclusion
In this paper, the discrete interface bond model is obtained by introducing the nonlinear shear bond strength criterion. e improved model solves the intrinsic problem of classical bonded particle model (BPM) in PFC. Based on this new model for simulating the high-strength gypsum material, numerical analysis of the ball-surface contact breakage experiment of the spherical sample with a diameter of 50 mm was carried out. is provides a basis for the analysis of ballto-surface contact breakage from a mesoscopic point of view and further verifies the validity and reliability of the model. e main conclusions are summarized as follows: (1) e tensile strength and shear strength of the specimens based on the nonlinear discrete interface bonded model are significantly higher than those of the classical bonded model, and this makes up for the deficiency that the shear strength envelope of the classical bonded model is linear. is can provide shear properties similar to those of actual materials. e mechanical properties of high-strength gypsum used in experiments can be accurately simulated by selecting model samples with reasonable mesoparameters.
(2) e macro-breakage morphology of numerical simulation reveals the generation of local breakage and the formation process of cone core from the mesoscale. e sample is finally broken into almost equal two parts. e numerical phenomenon is similar to the indoor test results.
(3) e force chain distribution results show that the contact in the conical area below the loading plate is subjected to a large compressive load. e angle between the force chain direction at the edge of the conical region and the vertical direction is approximately 45°. e maximum value of the tension chain is mainly distributed horizontally near the central longitudinal axis, and the distribution of the tension chain in the conical region is sparse. (4) Combined with the stress distribution and force chain distribution, the formation mechanism of cone core can be obtained. Under the action of the loading force, the conical area under the loading plate is under triaxial compression. When the contact force reaches the critical value, the shear stress on the inclined section reaches the shear strength of the material and causes shear failure. With the accumulation of local shear failure, the cone core is formed gradually.

Data Availability
All data used to support the findings of this study are included within the article.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this article.