Three-Dimensional Discrete Element Analysis of Crushing Characteristics of Calcareous Sand Particles

Department of Civil and Architecture Engineering, Jiangsu University of Science and Technology, Zhenjiang 212000, China Department of Civil Engineering, Suzhou University of Science and Technology, Suzhou 215009, China Department of Geotechnical Engineering, Nanjing Hydraulic Research Institute, Nanjing 210024, China College of Civil and Transportation Engineering, Hohai University, Nanjing, Jiangsu 210024, China


Introduction
Calcareous sand particles have the characteristics of low strength and are easy to break, which always occur obvious particle crushing [1,2]. If the particle crushing characteristic of calcareous sand particles is not considered in the design and construction process, the Island reef engineering construction will be directly affected [3]. Nowadays, a lot of studies on particle crushing of calcareous sand particles have been conducted by researchers.
Since the engineering accident happened in the northwest offshore drilling platform under developing in Australia in the 1980s, the particle crushing characteristics of calcareous sand have been widely studied by researchers around the world. For example, Wang et al. [4] conducted a series of drained and undrained triaxial compression test and estab-lished a constitutive model for calcareous sand considering particle crushing, which revealed the influence of particle crushing on deformation. Besides, Zhang and Luo [5] conducted a couple of triaxial tests which found that the particle crushing has important influence on the dilatancy and critical state of calcareous sand and proposed a modified dilatancy equation considering the effect of particle crushing. In addition, Wu et al. [6] conducted a series of drained and undrained triaxial compression tests under high confining pressure with different initial densities on the calcareous sand and determined multiple critical state lines with different initial densities. Furthermore, Wei et al. [7] performed multiple groups of ring shear tests under different initial conditions and analyzed the broken samples with SEM scanning and grading screening. Results found that the shear zone consists of highly crushed fine-grained calcareous sand, and the smoothness of particle surface increases with the increase of shear strain.
However, the production process of particle crushing cannot be controlled on real time by the current indoor test technology, and the contact mechanical properties between the inner particles cannot be tested [8,9]. In recent years, many researchers have achieved success in using discrete element to study particle breakage, which can supplement and expand the indoor test. For instance, based on bonded-particle model (BPM) and triaxial test, Wang et al. [10] found that the critical state of particle crushing is determined by the final particle gradation. In addition, through the simulation of one dimensional compression test on three fragment replacement modes, Zhou et al. [11] concluded that the particle crushing and the hardening law of particle strength may be the reason for the fractal of particle system. Besides, in the work of McDowell and de Bono [12], three fragment replacement modes are constructed with a certain amount of overlap and were placed on the initial boundaries of sample, and the micromechanical properties were studied on the basis of confining compression test. Furthermore, Ciantia et al. [13] adopted Apollonian filling method to construct multiple fragment replacement modes and conducted confining compression test. Results showed that the mechanical responses are significantly different under different replacement modes. Moreover, Cil and Buscarnera [14] considered two kinds of particle crushing modes in simulating the sand with different particle size based on Weibull distribution theory, and the study confirmed the interaction between the variability of particle strength and yield stress. What is more, Shi et al. [15] took the calcareous sand as simulated object and studied the influence of the principal stress and particle crushing on critical state of sand and compared the mechanical properties at critical state under the conditions of particle crushing or not. The study found that under different critical stress ratios, the mechanical properties of samples crushed or not have significant differences. Additionally, Tong et al. [16] simulated the particle crushing of geotechnical material by one dimensional compression specimen and proposed a prediction model for the particle crushing of geotechnical material. The results showed that the coordination number plays a leading role in the gradation evolution. Differently, through indoor test and threedimensional discrete element simulation, Lv et al. [17] studied the response of particle morphology and internal porosity of calcareous sand to particle mechanical properties and found that the particle strength decreases with the increase of porosity. Further, Kuang et al. [18] used BPM to analyze the single-particle strength of calcareous sand and studied its failure mode under different particle size. Therefore, nowadays, the discrete modeling studies on particle crushing of calcareous sand, whether BPM or FRM, can all feedback dynamically the real-time particle crushing, which is convenient for the in-depth study of the micromechanical mechanism of particle crushing.
In summary, on the basis of the indoor triaxial consolidation drainage test of calcareous sand and the numerical simulation platform of discrete element software PFC 3D , the meso model of calcareous sand particle crushing is estab-lished in this paper. In numerical simulation, FRM method is used to simulate the particle crushing, and the 14fragment replacement method satisfying Apollonian distribution is adopted. The evolution of macro and meso parameters describing the crushing characteristics of calcareous sand particles, such as porosity, effective porosity, effective coordination number, and sliding contact ratio, is analyzed, and further, the particle crushing mechanism of calcareous sand is studied.

Discrete Element Model for
Particle Crushing 2.1. Particle Crushing Criterion. de Bono and McDowell [19][20][21][22] systematically studied the particle crushing criterion, and the researchers indicated that the numerical tests were confirmed to the indoor test results and concluded that the octahedral stress criterion is suitable for sand particle crushing. In this paper, the particle crushing criterion proposed by de Bono and McDowell [19][20][21][22] is adopted as the crushing criterion, and the octahedral shear stress can be described as follows: where σ 1 , σ 2 , and σ 3 are the large, medium, and small principal stresses acted on any particle in the sample. The octahedral crushing criterion is proposed according to the homogenization thought, that is, in the numerical model, the inner stress distributes uniformly in any particles, and the inner stress can be expressed as Equation (2): where V P is the particle volume, N c is the contact number of particles, F ðcÞ j is the contact force at contact point ðcÞ along j direction, n ðcÞ i is the component along direction at contact point, and i, j are the Cartesian component.
In the simulation process in this paper, when the octahedral shear stress of one particle is greater than the set allowable shear stress in the sample, the particle is determined to be crushed. de Bono and McDowell [23] proposed the following relationship between the allowable shear stress and characteristic stress σ f obtained from single-particle crushing test: Jiang et al. [24] statistically analyzed the single-particle strength of calcareous sand and found that its crushing strength satisfies the Weibull distribution, which also has size effect. Based on Weibull distribution theory, the characteristic stress of each calcareous sand particle σ f can be expressed as follows: 2 Geofluids where P s ðdÞ, m, d 0 , and σ 0 represent the particle survival probability, Weibull modulus, reference particle size, and crushing strength corresponding to reference particle size.

Replacement Mode for Particle
Crushing. Based on the Apollonian filling method [25] in mathematics field, Ciantia et al. [13] constructed five particle crushing replacement modes and found that using the replacement method of 14 fragments can obtain the mechanical response similar to the replacement mode of 57 fragments. The Apollonian filling method requires that the subparticles in the parent particles in three-dimensional space are tangent to each other, and the subparticles are all tangent to the parent particles.
In order to insure the accuracy of numerical simulation, in this paper, the 14-fragment replacement mode satisfying the Apollonian distribution is adopted, as shown in Figure 1. The simulation process of this mode is as follows: for the target particles in the sample, when it satisfies the octahedral crushing criterion, delete the target particle and replace it with the 14 subparticles satisfying Apollonian distribution, and the distribution information of subparticles is listed in Table 1. In order to satisfy the principle of mass conservation and volume conservation, the subparticles are enlarged by 1.23 times by particle size expansion method, and in order to avoid local stress increase caused by subparticle expansion, the subroutine is compiled by FISH language, so that the subparticles can be amplified at low speed in local time steps.

Numerical Test Preparation
3.1. Preparation of Discrete Element Samples. In this paper, the numerical calculation platform is used to conduct the triaxial consolidated drained shear test. The same gradation to the indoor triaxial test is adopted (as shown in Figure 2), and the sample has a diameter of 39.1 mm and a height of 80 mm. The numerical test considers not the influence of particle shape and uses sphere basic units to represent calcareous sand, and the numerical samples are presented in Figure 3.
As the gradation of the sample has a wide range and many particle numbers, if the crushed subparticles satisfy the octahedral crushing criterion during calculation process, they will be further crushed into smaller subparticles, and this process is called the second crushing. Considering the calculation efficiency, the number of crushed subparticles should be controlled. Here, a crushing critical diameter is set d limit = 0:25d min , and if the particle diameter d is smaller than the critical particle diameter d limit , no more particle crushing will occur. The particle crushing samples in the numerical simulation are shown in Figure 4. It can be seen that the triaxial numerical sample has a diameter of 39.1 mm and a height of 80 mm. The particles are generated according to the gradation and the radius expansion method; first, the particle velocity is set to 0, and then, the uniform confining pressure is applied to the wall around the model through the servo mechanism. The interaction between particles adopts nonlinear Hertz-Mindlin contact model, and the particle crushing criterion parameters are set according to Reference [24]; the parameters in numerical model are listed in Table 2.
According to the test results in Reference [26], the particle meso parameters are determined by repeated trial   Note: R is the radius of parent particle. 3 Geofluids calculation of the numerical model. By comparison analysis, it can be seen that the mechanical responses of calcareous sand under triaxial shear condition can be well reflected in the numerical simulation results.

Scheme of Numerical Test.
During the triaxial test process, particle crushing occurs in both consolidation and shear stages. Based on previous studies, particle crushing mainly occurs in the shear stage; thus, in the present paper, the influence of consolidation on particle crushing is not considered. Besides, in order to compare the micromechanical difference with and without considering particle crushing, another group of simulation test considering not particle crushing is conducted. The numerical simulation scheme is shown in Table 3, with a total of 8 tests.

Relative Crushing Ratio Comparison of Numerical Test and Physical
Test. In order to quantify the degree of the particle crushing, Einav [27] modified the relative crushing rate B r proposed by Hardin [28], and the calculation formula is as follows: where d M is the maximum particle diameter; d m is the minimum particle diameter; F 0 ðdÞ, FðdÞ, and F u ðdÞ are the grading curve fitting function before loading, the current grading curve fitting function after loading, and the grading curve fitting function under limit state, respectively. The comparison of the relative crushing rate obtained by numerical simulation and indoor test is plotted in Figure 5. From this figure, it can be seen that the relative crushing rate in numerical simulation and indoor test both increase with the increase of confining pressure, and their difference becomes more significant. This is because that under high confining pressure condition, the inner stress in the sample is high, the inner particles are easy to crush, and second or multiple crush may occur. The final particle number of the sample under different working conditions is listed in Table 4. It can be seen that the particle shape also has an influence, the particles in the numerical samples are all spherical particles, while in the test, the shape of the calcareous sand is irregular. Therefore, in the test process, the obtained limit gradation is higher than that in the numerical simulation, which causes the B P * in the indoor test higher than that in the numerical simulation. Thus, in the same confining pressure, the simulated B P * is higher than that in the indoor test, and this regulation is in consistence with the References [10,15,29,30].

Influence of Particle Crushing on Mechanical
Coordination Number of Sample. The coordination number of particle [11,16,31,32] refers to the contact number of one particle to surrounding particles, which is an important parameter reflecting the inner characteristics of the sample.   4 Geofluids In order to remove the influence of suspended particles on the sample, the mechanical coordination number CN m defined by Thornton [33] is used to analyze the particle crushing on sample.
where N c is the total contact number of particles; N p , N p0 , and N p1 are the total number of particles, the number of particles with coordination number 0, and the number of particles with coordination number 1 in the sample. The mechanical properties of the calcareous sand before and after crushing can be well reflected by using the mechanical coordination number. As shown in Figure 6, under different confining pressure, the evolution of the coordination number of calcareous sand with and without considering particle crushing is similar. This is consistent with the research results of literature [34].
From Figure 6, it can be seen that under different confining pressure, the relationship between the mechanical coordination number and axial strain of sample considering particle crushing has a similar trend, which increases slowly first and then decreases sharply. Under different confining pressure, the mechanical coordinate number has significant difference, and the larger the confining pressure, the more the mechanical coordinate number, which indicates that the larger the confining pressure, the more the crushed particles in the sample, and the more contacts the crushed particles with neighboring particles, which is in consistence with the results in References [16,35,36].
Under the same confining pressure, the coordinate number of the sample considering not particle crushing is more than that considering the particle crushing, and with the increase of confining pressure, this trend is more significant. This is due to the fact that when particle crushing occurs, the internal structure of the sample becomes worse, and the greater the amount of particle crushing, the more obvious it is to the overall structure of the sample, resulting in the reduction of the stiffness and bearing capacity of the sample.

Influence of Particle Crushing on Porosity.
Porosity is an important index to describe the compactness degree of sample [37][38][39], and the evolution curves of porosity under different working conditions are plotted in Figure 7.
From Figure 7, it can be seen that the porosity of the sample decreases with the increase of confining pressure, Note: D is the sample diameter; H is the sample height; ρ s is the particle density; μ is the friction coefficient of particles; G is the shear modulus of particles; ν is the Poisson's ratio of particles; d limit is the ultimate particle size of particle breakage; d 0 is the reference particle size; σ 0 is the crushing strength corresponding to the reference particle size; m is the Weibull modulus. Note: in order to have convenient description, the samples are numbered, such as C100 and U100. C means the particle crushing is considered, U means the particle crushing is not considered, while 100 means the current confining pressure of the sample is 100 kPa. Other sample numbers can be inferred in the same way.    5 Geofluids and the evolution trend becomes gentle. This is because with the increase of the confining pressure, the sample becomes more dense on the whole, and as the sample shear shrinks first and then expands under low confining pressure, inducing the porosity first decreases and then increases slowly. Besides, as the calcareous sand is easy to crush, even under low confining pressure, significant particle crushing can also occur. Comparing with the evolution curve of porosity without considering particle crushing, it can be seen that particle crushing increases the compressibility of the sample, which also agrees with Reference [40].
During the test process, some small particles will dope in the gap between large particles but do not bear the external force. In order to describe the ratio of particles bearing force, the effective porosity n eff defined by Thornton [33] is used, i.e., the porosity of the sample after removing the suspended particles whose number of particles in contact which is less than 2.
where V t is the sample volume and V s and V f represent the volume of particles and the volume of suspended particles in the sample, respectively. From Figure 8, when the influence of suspended particles is removed and the axial strain is larger than 14%, the effective porosity increases slowly considering particle crushing. This indicates that at the later shear stage of sample, particles occur multistage crush, and the multistage crushed particles do not have good contact with the neighbor particles, and most of the filled particles are multistage crushed particles. As the larger the confining pressure, the more the crushed particles, the gap between particles is gradually filled by fine particles, and the particle crushing makes the sample more dense; thus, the porosity of samples C300 and C400 is lower than their effective porosity.

Influence of Particle Crushing on Sliding Contact Ratio.
During the shear process, the inner particles crush, inducing the redistribution of particles. In order to quantify this  6 Geofluids characteristic, the sliding contact ratio index (short for slip ratio) proposed by Bolton et al. [41] is used.
where N s is the sliding contact, i.e., the number of the contact number of the particle whose friction coefficient is greater than 99% (μ = 0:5) of the particle friction coefficient in the sample; N c indicates the number of contacts between particles in the sample. The evolution trend of the sliding contact ratio of each group of samples is plotted in Figures 9 and 10, from which it can be seen that with the increase of axial strain, the sliding ratio considering particle crushing increases first and then decreases and becomes stale when axial strain is larger than 13%. This is because at the early and middle stages of test (axial stain ≤ 8%), the increasing vertical stress of the sample leads to a large number of particle crushing, rotation,  Axial strain, 1 (%) U100 U200 U300 U400 Figure 10: Evolution of sliding contact ratio without considering particle breakage. 7 Geofluids and rolling and further leads to the increase of the number of particles and contact. The particle breakage tends to be stable in the middle and later stages of the test (axial strain 8%), and the contact force distribution of the sample is more uniform.
For the samples without considering particle crushing, with the increase of axial strain and confining pressure, the sliding ratio changes significantly, which is due to the fact that during the loading process, the gap between large particles is not filled with medium and small particles after crushing, and there are more particles in the sliding state of the sample, which is not conducive to the redistribution of particle position. Therefore, particle breakage will make the properties of calcareous sand more uniform, which is in consistence with the results in Reference [42].

Conclusions
Based on the Apollonian filling algorithm and indoor test and under the premise of satisfying the conservation of particle mass and no particle overlap, a three-dimensional discrete element model considering particle breakage was established, and the influence of particle crushing on the macro and meso mechanical properties of calcareous sand was analyzed. The following conclusions can be drawn: (1) The numerical simulation reproduces the indoor test results and is in consistence with the existed research findings. Thus, the present numerical model can provide a new method for analyzing particle crushing of calcareous sand (2) With the increase of confining pressure, the relative crushing ratio of sample increases accordingly, and with the increase of confining pressure, the difference of the crushing ratio between the numerical simulation and indoor test becomes larger. Thus, during the simulation process, using different particle shape fragment replacement model is a problem to be solved in fine numerical simulation (3) The mechanical coordination number under different confining pressures increases first and then decreases with the shear process, and the lower the confining pressure, the more obvious the trend is.
In the shear process, the dilatancy of the specimen under low confining pressure is one of the reasons for the instability of the mechanical coordination number (4) There are significant differences in micro mechanical properties between the samples with and without particle crushing. The particle crushing makes the material properties of the samples more uniform and stable

Data Availability
The data used to support the findings of this study are included in the article.

Conflicts of Interest
There is no conflict of interest in this study.