Micro-Mechanism of Spherical Gypsum Particle Breakage under Ball–Plane Contact Condition

: Coarse-grained soils are used extensively in engineering applications. The breakage of coarse-grained soil particles may have a great e ﬀ ect on their mechanical characteristics. It is important to fully understand the phenomenon of particle breakage and comprehend its e ﬀ ect on engineering properties. The aim of this study was to investigate the process and mechanism of spherical particle breakage under ball–plane contact conditions. Particle contact tests and corresponding simulations based on the discrete element method were performed. The mechanical properties and breaking morphologies of gypsum balls, as well as the signiﬁcant feature of the existence of a cone core under the contact point, were obtained by the experiments. To enable particle crushing in a numerical simulation, noncrushable elementary particles were bonded together to represent the specimen. The numerical model, which was validated by the unconﬁned compression test and splitting test, was well ﬁtted with the experiment by applying ﬂat-joint contact. More importantly, the combination of the simulation and experiment demonstrated the role that the cone core plays during particle breakage and revealed the mechanism of the formation of the cone core and its e ﬀ ect on particle breakage.


Introduction
Coarse-grained soils exhibit good compaction performance, high strength, strong permeability, and high bearing capacity. It is generally acknowledged that the mechanical behavior of coarse-grained soils is very complex. The crushing mechanism of coarse-grained soils and their behavior under high stress has garnered the attention of scholars. It is important to thoroughly understand the phenomenon of particle breakage and comprehend its effect on engineering properties [1]. A large number of experiments and numerical studies have been carried out on the crushing and mechanical behavior of coarse-grained soil particles. However, the mechanism of particle breakage needs further study.
At present, the main methods used to experimentally study particle breakage include the triaxial test [2][3][4][5][6][7][8][9][10], the shear test [11][12][13][14][15], the uniaxial compression test [8,9,16], the particle contact test [17,18] and the centrifuge modeling test [19]. However, previous studies have primarily focused on the behavior of soil aggregates. The physical and mechanical properties of particle breakage have been examined from the perspective of the average stress of the aggregate, but the mechanical behavior of a single particle has rarely been analyzed. The failure processes and fracture patterns of single particle breakage have been analyzed in recent years. However, it is still difficult to observe the evolution of crushing or reveal the basic mechanism of particle breakage through a conventional experiment.
The results obtained by experiment and simulation demonstrate the formation of the cone core and reveal its effect on particle breakage.

Experimental Study
Considering that it is difficult to carry out particle-contact tests under various conditions at one time, particle contact was simplified to a ball-plane contact in this study. Moreover, the ball-plane contact was simulated by a spherical particle-rigid wall contact. The breakage behavior of spherical specimens under ball-plane contact conditions was studied by repeated experiments.

Material and Methods
Due to its advantage in sample preparation, high-strength gypsum was selected as the test material. In the experiment, spherical particles of high-strength gypsum with a diameter of 50 mm were selected. If the particle size is too large, the exothermic hardening process of gypsum may induce temperature cracking inside the specimen. However, if the particle size is too small, there could be difficulties in observing the morphology after breakage. Therefore, only 50 mm particles were produced.
Unconfined compression tests and splitting tests were performed to measure the basic mechanical parameters for high-strength gypsum material. The specimens for compression tests were 50 mm in diameter and 100 mm in height (see Figure 1a). The specimens for splitting tests were 50 mm in diameter and 50 mm in height (see Figure 1b). The basic mechanical parameters acquired by rock tests for high-strength gypsum material are shown in Table 1. The results are the means of 5 repeated tests.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 3 of 16 Considering that it is difficult to carry out particle-contact tests under various conditions at one time, particle contact was simplified to a ball-plane contact in this study. Moreover, the ball-plane contact was simulated by a spherical particle-rigid wall contact. The breakage behavior of spherical specimens under ball-plane contact conditions was studied by repeated experiments.

Material and Methods
Due to its advantage in sample preparation, high-strength gypsum was selected as the test material. In the experiment, spherical particles of high-strength gypsum with a diameter of 50 mm were selected. If the particle size is too large, the exothermic hardening process of gypsum may induce temperature cracking inside the specimen. However, if the particle size is too small, there could be difficulties in observing the morphology after breakage. Therefore, only 50 mm particles were produced.
Unconfined compression tests and splitting tests were performed to measure the basic mechanical parameters for high-strength gypsum material. The specimens for compression tests were 50 mm in diameter and 100 mm in height (see Figure 1a). The specimens for splitting tests were 50 mm in diameter and 50 mm in height (see Figure 1b). The basic mechanical parameters acquired by rock tests for high-strength gypsum material are shown in Table 1. The results are the means of 5 repeated tests.  The experiments are carried out by employing a rock rheological testing system (as shown in Figure 2), which can provide stable displacement-controlled and servo-controlled loading conditions. The spherical joint support is specifically customized to prevent eccentric compression during the test (see Figure 2b). The curvature radius of the support was much larger than the radius of the specimens. The precision of measurement and control for stress and deformation were 5 N and 0.001 mm, respectively. The frequency of data acquisition was 1 Hz. The ball-plane contact condition of the test is shown in Figure 3. The displacement and contact force of the loading equipment were recorded by the testing system. Additional extensometers placed in front and at the back of the specimen were also used to measure the displacement of the supports.  The experiments are carried out by employing a rock rheological testing system (as shown in Figure 2), which can provide stable displacement-controlled and servo-controlled loading conditions. The spherical joint support is specifically customized to prevent eccentric compression during the test (see Figure 2b). The curvature radius of the support was much larger than the radius of the specimens. The precision of measurement and control for stress and deformation were 5 N and 0.001 mm, respectively. The frequency of data acquisition was 1 Hz. The ball-plane contact condition of the test is shown in Figure 3. The displacement and contact force of the loading equipment were recorded by the testing system. Additional extensometers placed in front and at the back of the specimen were also used to measure the displacement of the supports.  The loading progress was displacement controlled, and the contact force-displacement was measured during the test. The loading rate was set to 0.002 mm/s. Though a ball-plane contact test does not evaluate strain rate, the loading rate divided by the diameter of the specimen belongs to a quasi-static procedure [42]. Before the tests started, the loading plates and fixture were lubricated with oil to eliminate the effect of friction. The measuring system was reset before every test. Approximately 200 N were applied to ensure proper seating between the loading fixtures and the specimen.

Mechanical Characteristics
The force displacement curves for the five tested specimens are shown in Figure 4. The displacement is the average of the data obtained by the sensors placed around the specimen. Force was obtained from the load cell of the loading system.
As shown in the figure, the curves are almost linear. The initial slightly concave nature may indicate local breakage. Initially, local breakage occurred at the contact point due to stress concentration, and the equivalent stiffness was relatively small, at approximately 2.871 kN/mm (see A in Figure 4). The slope of the force-displacement curve at this stage is almost unchanged, and the equivalent stiffness was approximately 7.953 kN/mm. The last stage was particle breakage. By the time the contact force exceeded its strength, the gypsum ball suddenly crushed into blocks.
The average maximum loading force for the ball-plane contact of gypsum particles was approximately 6.843 kN. Specimen 1 was partially crushed just before the overall breakage. Specimen 3 was crushed at a relatively small contact force because of manufacturing defects. A bubble inside Specimen 3 led to this result, and it may have been formed during the casting process of gypsum.  The loading progress was displacement controlled, and the contact force-displacement was measured during the test. The loading rate was set to 0.002 mm/s. Though a ball-plane contact test does not evaluate strain rate, the loading rate divided by the diameter of the specimen belongs to a quasi-static procedure [42]. Before the tests started, the loading plates and fixture were lubricated with oil to eliminate the effect of friction. The measuring system was reset before every test. Approximately 200 N were applied to ensure proper seating between the loading fixtures and the specimen.

Mechanical Characteristics
The force displacement curves for the five tested specimens are shown in Figure 4. The displacement is the average of the data obtained by the sensors placed around the specimen. Force was obtained from the load cell of the loading system.
As shown in the figure, the curves are almost linear. The initial slightly concave nature may indicate local breakage. Initially, local breakage occurred at the contact point due to stress concentration, and the equivalent stiffness was relatively small, at approximately 2.871 kN/mm (see A in Figure 4). The slope of the force-displacement curve at this stage is almost unchanged, and the equivalent stiffness was approximately 7.953 kN/mm. The last stage was particle breakage. By the time the contact force exceeded its strength, the gypsum ball suddenly crushed into blocks.
The average maximum loading force for the ball-plane contact of gypsum particles was approximately 6.843 kN. Specimen 1 was partially crushed just before the overall breakage. Specimen 3 was crushed at a relatively small contact force because of manufacturing defects. A bubble inside Specimen 3 led to this result, and it may have been formed during the casting process of gypsum. The loading progress was displacement controlled, and the contact force-displacement was measured during the test. The loading rate was set to 0.002 mm/s. Though a ball-plane contact test does not evaluate strain rate, the loading rate divided by the diameter of the specimen belongs to a quasi-static procedure [42]. Before the tests started, the loading plates and fixture were lubricated with oil to eliminate the effect of friction. The measuring system was reset before every test. Approximately 200 N were applied to ensure proper seating between the loading fixtures and the specimen.

Mechanical Characteristics
The force displacement curves for the five tested specimens are shown in Figure 4. The displacement is the average of the data obtained by the sensors placed around the specimen. Force was obtained from the load cell of the loading system.
As shown in the figure, the curves are almost linear. The initial slightly concave nature may indicate local breakage. Initially, local breakage occurred at the contact point due to stress concentration, and the equivalent stiffness was relatively small, at approximately 2.871 kN/mm (see A in Figure 4). The slope of the force-displacement curve at this stage is almost unchanged, and the equivalent stiffness was approximately 7.953 kN/mm. The last stage was particle breakage. By the time the contact force exceeded its strength, the gypsum ball suddenly crushed into blocks.
The average maximum loading force for the ball-plane contact of gypsum particles was approximately 6.843 kN. Specimen 1 was partially crushed just before the overall breakage. Specimen 3 was crushed at a relatively small contact force because of manufacturing defects. A bubble inside Specimen 3 led to this result, and it may have been formed during the casting process of gypsum.

Breakage Behavior
The particle contact test process and its corresponding force-displacement curve are shown in Figure 5. The experimental results showed that all the specimens broke in a similar manner. We can see the typical features of gypsum particles after breakage in Figure 6.
At the start of the test, stress concentration at the contact between the specimen and the loading platen caused local breakage. The cracks first appeared in the tests when the contact force reached about 95% of its peak value. As the loading process continued, the crack propagated rapidly just before the overall breakage. The final breakage was a brittle fracture that the contact force instantaneously dropped to zero. By observing the morphologies of specimen after breakage, a cone core could be found under the contact area. The contact region is shown as Figure 6a, the cone area on the fracture surface is shown in Figure 6b, and the detail of the enlarged scale of the cone core can be seen in Figure 6c. The cone core was approximately 11 mm in diameter and 12 mm in depth. It is obvious that the cone core

Breakage Behavior
The particle contact test process and its corresponding force-displacement curve are shown in Figure 5. The experimental results showed that all the specimens broke in a similar manner. We can see the typical features of gypsum particles after breakage in Figure 6.
At the start of the test, stress concentration at the contact between the specimen and the loading platen caused local breakage. The cracks first appeared in the tests when the contact force reached about 95% of its peak value. As the loading process continued, the crack propagated rapidly just before the overall breakage. The final breakage was a brittle fracture that the contact force instantaneously dropped to zero.

Breakage Behavior
The particle contact test process and its corresponding force-displacement curve are shown in Figure 5. The experimental results showed that all the specimens broke in a similar manner. We can see the typical features of gypsum particles after breakage in Figure 6.
At the start of the test, stress concentration at the contact between the specimen and the loading platen caused local breakage. The cracks first appeared in the tests when the contact force reached about 95% of its peak value. As the loading process continued, the crack propagated rapidly just before the overall breakage. The final breakage was a brittle fracture that the contact force instantaneously dropped to zero. By observing the morphologies of specimen after breakage, a cone core could be found under the contact area. The contact region is shown as Figure 6a, the cone area on the fracture surface is shown in Figure 6b, and the detail of the enlarged scale of the cone core can be seen in Figure 6c. The cone core was approximately 11 mm in diameter and 12 mm in depth. It is obvious that the cone core By observing the morphologies of specimen after breakage, a cone core could be found under the contact area. The contact region is shown as Figure 6a, the cone area on the fracture surface is shown in Figure 6b, and the detail of the enlarged scale of the cone core can be seen in Figure 6c. The cone core was approximately 11 mm in diameter and 12 mm in depth. It is obvious that the cone core penetrated into the ball particle during the test, which may have contributed to the final breakage of the particle. The formation process of the cone core could not be revealed by conventional experiments.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 6 of 16 penetrated into the ball particle during the test, which may have contributed to the final breakage of the particle. The formation process of the cone core could not be revealed by conventional experiments.

Numerical Simulation
The particle-contact test for gypsum was simulated by the three-dimensional particle flow code (PFC 3D ) (Itasca Consulting Group, Inc. Minneapolis, MN, USA), which was developed by Itasca and a widely used the DEM software program. The material-modeling support of the PFC was used for generating the specimen and for performing standard rock tests.
The numerical model in the PFC is composed of rigid balls that interact at their contacts. The PFC model was originally developed to simulate the micromechanical behavior of noncohesive media such as soils and sands [33]. In addition, the built-in bonded-particle model (BPM) for rock applies cohesive bonds between elementary balls to simulate the behaviors of solid rocks [35], in which crack formation is simulated through the breaking of bonds, while fracture propagation is obtained by the coalescence of multiple bond breakages. In this study, we bound uncrushable elementary balls together to represent specimens that can be crushed when the external force exceeds their strength.

Flat-Joint Model
A typical flat-joint contact can be represented by the model shown in Figure 7. A flat-jointed material consists of rigid balls joined by flat-joint contacts. The effective surface of each ball is defined by the notional surfaces of its pieces, which interact at each flat-joint contact with the notional surface of the contact piece. Balls in a flat-jointed material are depicted as spherical cores and several skirted faces, which can provide resistance to rotation even after bond breaking. The skirted faces are deleted when the relative displacement at a flat-joint contact becomes greater than the flat-joint diameter. If the two balls are in contact again, the behavior will be that of between spherical surfaces. The interface between notional surfaces is a line (in 2D) or a disk (in 3D) and discretized into a number of elements in both circumferential and radial directions. Since the destruction of each element is relatively independent, partial damage of a contact can be achieved.

Numerical Simulation
The particle-contact test for gypsum was simulated by the three-dimensional particle flow code (PFC 3D ) (Itasca Consulting Group, Inc. Minneapolis, MN, USA), which was developed by Itasca and a widely used the DEM software program. The material-modeling support of the PFC was used for generating the specimen and for performing standard rock tests.
The numerical model in the PFC is composed of rigid balls that interact at their contacts. The PFC model was originally developed to simulate the micromechanical behavior of noncohesive media such as soils and sands [33]. In addition, the built-in bonded-particle model (BPM) for rock applies cohesive bonds between elementary balls to simulate the behaviors of solid rocks [35], in which crack formation is simulated through the breaking of bonds, while fracture propagation is obtained by the coalescence of multiple bond breakages. In this study, we bound uncrushable elementary balls together to represent specimens that can be crushed when the external force exceeds their strength.

Flat-Joint Model
A typical flat-joint contact can be represented by the model shown in Figure 7. A flat-jointed material consists of rigid balls joined by flat-joint contacts. The effective surface of each ball is defined by the notional surfaces of its pieces, which interact at each flat-joint contact with the notional surface of the contact piece. Balls in a flat-jointed material are depicted as spherical cores and several skirted faces, which can provide resistance to rotation even after bond breaking. The skirted faces are deleted when the relative displacement at a flat-joint contact becomes greater than the flat-joint diameter. If the two balls are in contact again, the behavior will be that of between spherical surfaces. The interface between notional surfaces is a line (in 2D) or a disk (in 3D) and discretized into a number of elements in both circumferential and radial directions. Since the destruction of each element is relatively independent, partial damage of a contact can be achieved.
Each element can be either bonded or unbonded, and the breakage of each element contributes partial damage to the interface. The behavior of a bonded element is linearly elastic, and the bond breaks when the strength limit is exceeded. The behavior of an unbonded element is linearly elastic and frictional, and the shear strength meets the Coulomb criterion. A tensile load is carried only in a bonded region. The force-displacement response of the flat-joint interface includes evolving from a fully bonded state to a frictional state. Appl. Sci. 2019, 9,   Each element can be either bonded or unbonded, and the breakage of each element contributes partial damage to the interface. The behavior of a bonded element is linearly elastic, and the bond breaks when the strength limit is exceeded. The behavior of an unbonded element is linearly elastic and frictional, and the shear strength meets the Coulomb criterion. A tensile load is carried only in a bonded region. The force-displacement response of the flat-joint interface includes evolving from a fully bonded state to a frictional state.

Model Calibration and Validation
To simulate the macroscopic phenomena of the experimental test in numerical modeling, repeated calculation and adjustment were required for model calibration. The numerical model was constructed and calibrated against the experimental results. In this paper, parameter adjusting methods were referred to the PFC official handbook. Studies have shown that the size of elementary particles has a certain effect on the mechanical properties of the numerical model [43,44]. When the elementary particle size is too large, the macroscopic mechanical properties of the model are unstable. The macroscopic mechanical properties of the model are not affected when the size of the elementary particle is sufficiently small. When the elementary particle size is too small, the number of particles in a numerical model will greatly increase to a point where the computational complexity will be unacceptable. We set the minimum diameter of elementary particles to 1.0 mm and the ratio of maximum to minimum size to 1.6. The ratio of model diameter to median particle size ( / ) is 38.5, and the effect on the simulation is very small, according to Ding [45].
The unconfined compression test (UCT) and splitting test (ST) were used to calibrate and validate the microparameters of the model material. For comparison with the experiment, the numerical model was set to be the same size as the experimental specimen (see Figure 8). Each contact was initially bonded, gapped or slit (see Figure 9) in a flat-jointed material, where is the initial gap between the notional surfaces. The material microstructure parameters are the

Model Calibration and Validation
To simulate the macroscopic phenomena of the experimental test in numerical modeling, repeated calculation and adjustment were required for model calibration. The numerical model was constructed and calibrated against the experimental results. In this paper, parameter adjusting methods were referred to the PFC official handbook. Studies have shown that the size of elementary particles has a certain effect on the mechanical properties of the numerical model [43,44]. When the elementary particle size is too large, the macroscopic mechanical properties of the model are unstable. The macroscopic mechanical properties of the model are not affected when the size of the elementary particle is sufficiently small. When the elementary particle size is too small, the number of particles in a numerical model will greatly increase to a point where the computational complexity will be unacceptable. We set the minimum diameter of elementary particles to 1.0 mm and the ratio of maximum to minimum size to 1.6. The ratio of model diameter to median particle size (L/d) is 38.5, and the effect on the simulation is very small, according to Ding [45].
The unconfined compression test (UCT) and splitting test (ST) were used to calibrate and validate the microparameters of the model material. For comparison with the experiment, the numerical model was set to be the same size as the experimental specimen (see Figure 8). Each element can be either bonded or unbonded, and the breakage of each element contributes partial damage to the interface. The behavior of a bonded element is linearly elastic, and the bond breaks when the strength limit is exceeded. The behavior of an unbonded element is linearly elastic and frictional, and the shear strength meets the Coulomb criterion. A tensile load is carried only in a bonded region. The force-displacement response of the flat-joint interface includes evolving from a fully bonded state to a frictional state.

Model Calibration and Validation
To simulate the macroscopic phenomena of the experimental test in numerical modeling, repeated calculation and adjustment were required for model calibration. The numerical model was constructed and calibrated against the experimental results. In this paper, parameter adjusting methods were referred to the PFC official handbook. Studies have shown that the size of elementary particles has a certain effect on the mechanical properties of the numerical model [43,44]. When the elementary particle size is too large, the macroscopic mechanical properties of the model are unstable. The macroscopic mechanical properties of the model are not affected when the size of the elementary particle is sufficiently small. When the elementary particle size is too small, the number of particles in a numerical model will greatly increase to a point where the computational complexity will be unacceptable. We set the minimum diameter of elementary particles to 1.0 mm and the ratio of maximum to minimum size to 1.6. The ratio of model diameter to median particle size ( / ) is 38.5, and the effect on the simulation is very small, according to Ding [45].
The unconfined compression test (UCT) and splitting test (ST) were used to calibrate and validate the microparameters of the model material. For comparison with the experiment, the numerical model was set to be the same size as the experimental specimen (see Figure 8). Each contact was initially bonded, gapped or slit (see Figure 9) in a flat-jointed material, where is the initial gap between the notional surfaces. The material microstructure parameters are the Each contact was initially bonded, gapped or slit (see Figure 9) in a flat-jointed material, where g 0 is the initial gap between the notional surfaces. The material microstructure parameters are the fraction of bonded (φ B ) and gapped (φ G ) contacts. We assumed that there is no initial crack in an artificial gypsum specimen; thus, all the contacts were set to Type B. The interfaces were discretized into three elements in the circumferential direction (N α = 3), and there was one element in the radial direction (N r = 1). The radius multiplier was set to be fixed (C λ = 0) so that all flat-joint contact radius multipliers were set equal to the specific value (λ v = 1.0). The flat-joint contact model was set to be installed at all ball-ball contacts with a gap less than or equal to 0.0002 m (g i = 2.0 × 10 −4 m). The porosity-considered bulk density was set to 2650 kg/m 3 . The above parameters were fixed during the calibration process. fraction of bonded ( ) and gapped ( ) contacts. We assumed that there is no initial crack in an artificial gypsum specimen; thus, all the contacts were set to Type B. The interfaces were discretized into three elements in the circumferential direction (  3), and there was one element in the radial direction ( 1). The radius multiplier was set to be fixed ( 0) so that all flat-joint contact radius multipliers were set equal to the specific value ( 1.0). The flat-joint contact model was set to be installed at all ball-ball contacts with a gap less than or equal to 0.0002 m ( 2.0 × 10 m). The porosity-considered bulk density was set to 2650 kg/m . The above parameters were fixed during the calibration process. The calibration process involves the trial-and-error adjustment of the following micromechanical parameters: effective modulus * (flat joint) and * (linear), normal-to-shear stiffness ratio * (flat joint) and * (linear), friction coefficient (flat joint) and (linear), tensile strength , cohesion , and friction angle . All the tested values of these micromechanical parameters are listed in Table 2, and the calibration work was conducted using orthogonal experimental design. Note that in this study, * = * , * = * , and = .
The macroproperties used for model calibration include unconfined compressive strength, tensile strength, elastic modulus, and Poisson's ratio. The model sensitivity to key micromechanical parameters was analyzed with a multifactor Analysis of Variance (ANOVA) at 5% significance level ( 0.05). The F values and associated p-values are shown in Figure 10. We can accept the hypothesis that the property is sensitive to the factor if the corresponding p-value is less than ; otherwise, the property and the factor are relatively independent variates. The calibration process involves the trial-and-error adjustment of the following micromechanical parameters: effective modulus E * (flat joint) and E * lnm (linear), normal-to-shear stiffness ratio κ * (flat joint) and κ * lnm (linear), friction coefficient µ (flat joint) and µ lnm (linear), tensile strength σ c , cohesion c, and friction angle φ. All the tested values of these micromechanical parameters are listed in Table 2, and the calibration work was conducted using orthogonal experimental design. Note that in this study, E * = E * lnm , κ * = κ * lnm , and µ = µ lnm .  The normal and shear stiffnesses were set based on a specified deformability by defining the effective modulus (E * ) and normal to shear stiffness ratio (κ * ): The macroproperties used for model calibration include unconfined compressive strength, tensile strength, elastic modulus, and Poisson's ratio. The model sensitivity to key micromechanical parameters was analyzed with a multifactor Analysis of Variance (ANOVA) at 5% significance level (α = 0.05).
The F values and associated p-values are shown in Figure 10. We can accept the hypothesis that the property is sensitive to the factor if the corresponding p-value is less than α; otherwise, the property and the factor are relatively independent variates.
The following correlations are evident in Figure 10: 1.
The unconfined compressive strength is sensitive to σ c , as well as c and φ; 2.
The tensile strength is highly sensitive to σ c and can be affected by κ * ; 3.
The elastic modulus is highly sensitive to E * and κ * ; 4.
Poisson's ratio is sensitive to κ * and will be affected by c and φ.
During the model calibration process, some tests were applied with different random seeds to check the model's reliability. The results showed that the calculated coefficients of variation for the simulated macroproperties were less than 1%.
The above observations indicate certain correlations between the micromechanical parameters and macroproperties. However, this paper does not intend to establish a quantitative prediction model of the macroproperties of flat-jointed material. The parametric study was primarily conducted to optimize the calibration process. A set of best-fit parameters were obtained based on the parametric study, as listed in Table 3. Note that in this study, E * = E * lnm , κ * = κ * lnm , and µ = µ lnm .  The following correlations are evident in Figure 10: 1. The unconfined compressive strength is sensitive to , as well as and ; 2. The tensile strength is highly sensitive to and can be affected by * ; 3. The elastic modulus is highly sensitive to * and * ; 4. Poisson's ratio is sensitive to * and will be affected by and .
During the model calibration process, some tests were applied with different random seeds to check the model's reliability. The results showed that the calculated coefficients of variation for the simulated macroproperties were less than 1%.
The above observations indicate certain correlations between the micromechanical parameters and macroproperties. However, this paper does not intend to establish a quantitative prediction model of the macroproperties of flat-jointed material. The parametric study was primarily conducted to optimize the calibration process. A set of best-fit parameters were obtained based on the parametric study, as listed in Table 3. Note that in this study, * = * , * = * , and = .

Parameters Value
The test results were compared with those of the experiments. The force-displacement curves are shown in Figure 11. The result for the experiment was the average curve of five repeated tests.
The displacement during the initial stage of the splitting test was relatively large. A possible reason for this was that the contact area was relatively small and the steel bar used for loading may not have been absolutely flat; it was difficult to ensure complete contact between a specimen and the loading equipment; thus, the initial stiffness was relatively small. However, the slope and peak value of the numerical model corresponded with those of the experiment very well.
The test results were compared with those of the experiments. The force-displacement curves are shown in Figure 11. The result for the experiment was the average curve of five repeated tests. The displacement during the initial stage of the splitting test was relatively large. A possible reason for this was that the contact area was relatively small and the steel bar used for loading may not have been absolutely flat; it was difficult to ensure complete contact between a specimen and the loading equipment; thus, the initial stiffness was relatively small. However, the slope and peak value of the numerical model corresponded with those of the experiment very well.
The results show that the numerical model was well fitted with the mechanical properties of high-strength gypsum we used and could reflect its basic mechanical characteristics. The microparameters were well calibrated and were suitable for the study of particle breakage under ball-plane contact conditions. The physical characteristics for numerical material are shown in Table 4  The shape and size of the numerical model were consistent with the spherical gypsum particles in the experimental test. The diameter of the model was 50 mm. The minimum diameter of elementary particles was 1.0 mm, and the ratio of maximum to minimum size was 1.6. The ratio of model diameter to median particle size ( / ) was 38.5, and the effect on the simulation is generally very small, according to Ding [45].
The original model is shown in Figure 12. We used the wall unit in the PFC to simulate the contact surface and the loading plate. The loading speed of the numerical simulation was set in reference to that of the experiment. Note that the contact condition at the bottom was different from that of the experiment. According to Saint-Venant's Principle, the difference between the effects of different contact conditions at the bottom on the upper part of the specimen is very small. In this paper, only the morphology of the upper part is discussed to study the particle breakage under ballplane contact conditions. The results show that the numerical model was well fitted with the mechanical properties of high-strength gypsum we used and could reflect its basic mechanical characteristics. The microparameters were well calibrated and were suitable for the study of particle breakage under ball-plane contact conditions. The physical characteristics for numerical material are shown in Table 4  The shape and size of the numerical model were consistent with the spherical gypsum particles in the experimental test. The diameter of the model was 50 mm. The minimum diameter of elementary particles was 1.0 mm, and the ratio of maximum to minimum size was 1.6. The ratio of model diameter to median particle size (L/d) was 38.5, and the effect on the simulation is generally very small, according to Ding [45].
The original model is shown in Figure 12. We used the wall unit in the PFC to simulate the contact surface and the loading plate. The loading speed of the numerical simulation was set in reference to that of the experiment. Note that the contact condition at the bottom was different from that of the experiment. According to Saint-Venant's Principle, the difference between the effects of different contact conditions at the bottom on the upper part of the specimen is very small. In this paper, only the morphology of the upper part is discussed to study the particle breakage under ball-plane contact conditions.
We also used the built-in FISH functions in the PFC to enable fragment visualization by painting broken contacts with different colors, which made it possible to observe the particle breakage behavior. Appl. Sci. 2019, 9, x FOR PEER REVIEW 11 of 16

Figure 12. Original numerical model
We also used the built-in FISH functions in the PFC to enable fragment visualization by painting broken contacts with different colors, which made it possible to observe the particle breakage behavior.

Results
In the numerical simulation, the wall units worked as loading planes and measured the displacement as well as the contact force. The main purpose of this paper is not to validate the statistical rules of maximum contact force, which has been the subject of many previous studies [22][23][24]26,27]. The DEM simulation was intended to determine the general mechanism of contact breakage under ball-plane conditions from a microscopic perspective.

Mechanical Characteristics
The force-displacement curve recorded by numerical simulation was compared with that of the experiments, as shown in Figure 13. The numerical model fit the force-displacement curve very well in the particle contact test. The maximum ball-plane contact force of the numerical model was 6.989 kN, and the equivalent stiffness was approximately 8.160 kN/mm. It can also be seen from the curve that the local breakage of the numerical result was not evident at the initial stage. The reason could be that the model was an aggregate of a finite number of rigid elementary particles, which limited the degree of fragmentation at a relatively small displacement. However, the advantage of this model on the simulation of particle contact was appreciable when it was compared with previous studies using a parallel-bond model [18], the result of which had an unrealistically large displacement.

Results
In the numerical simulation, the wall units worked as loading planes and measured the displacement as well as the contact force. The main purpose of this paper is not to validate the statistical rules of maximum contact force, which has been the subject of many previous studies [22][23][24]26,27]. The DEM simulation was intended to determine the general mechanism of contact breakage under ball-plane conditions from a microscopic perspective.

Mechanical Characteristics
The force-displacement curve recorded by numerical simulation was compared with that of the experiments, as shown in Figure 13. The numerical model fit the force-displacement curve very well in the particle contact test. The maximum ball-plane contact force of the numerical model was 6.989 kN, and the equivalent stiffness was approximately 8.160 kN/mm. It can also be seen from the curve that the local breakage of the numerical result was not evident at the initial stage. The reason could be that the model was an aggregate of a finite number of rigid elementary particles, which limited the degree of fragmentation at a relatively small displacement. However, the advantage of this model on the simulation of particle contact was appreciable when it was compared with previous studies using a parallel-bond model [18], the result of which had an unrealistically large displacement. We also used the built-in FISH functions in the PFC to enable fragment visualization by painting broken contacts with different colors, which made it possible to observe the particle breakage behavior.

Results
In the numerical simulation, the wall units worked as loading planes and measured the displacement as well as the contact force. The main purpose of this paper is not to validate the statistical rules of maximum contact force, which has been the subject of many previous studies [22][23][24]26,27]. The DEM simulation was intended to determine the general mechanism of contact breakage under ball-plane conditions from a microscopic perspective.

Mechanical Characteristics
The force-displacement curve recorded by numerical simulation was compared with that of the experiments, as shown in Figure 13. The numerical model fit the force-displacement curve very well in the particle contact test. The maximum ball-plane contact force of the numerical model was 6.989 kN, and the equivalent stiffness was approximately 8.160 kN/mm. It can also be seen from the curve that the local breakage of the numerical result was not evident at the initial stage. The reason could be that the model was an aggregate of a finite number of rigid elementary particles, which limited the degree of fragmentation at a relatively small displacement. However, the advantage of this model on the simulation of particle contact was appreciable when it was compared with previous studies using a parallel-bond model [18], the result of which had an unrealistically large displacement.

Breakage Behavior
In the numerical test, the breakage phenomenon of the specimen was consistent with the experimental result. First, local breakage was produced, and a small platform was formed at the contact point. Then, the stage for crack nucleation and crack propagation followed. Finally, the contact force exceeded its strength, and the specimen broke into two pieces.
By means of numerical simulation, we could observe the law of morphological changes inside the specimen as well as the specific form of failure. The top view of the overall breakage shows the destruction form and the cracks in the surface layer of the specimen (see Figure 14a). The elementary particles and failed bonds are shown in the figure. The tensile and shear failure of bonds are indicated by arrows and texts (see Figure 14b and Figure 16c).

Breakage Behavior
In the numerical test, the breakage phenomenon of the specimen was consistent with the experimental result. First, local breakage was produced, and a small platform was formed at the contact point. Then, the stage for crack nucleation and crack propagation followed. Finally, the contact force exceeded its strength, and the specimen broke into two pieces.
By means of numerical simulation, we could observe the law of morphological changes inside the specimen as well as the specific form of failure. The top view of the overall breakage shows the destruction form and the cracks in the surface layer of the specimen (see Figure 14a). The elementary particles and failed bonds are shown in the figure. The tensile and shear failure of bonds are indicated by arrows and texts (see Figure 14b and 16c). The result indicates that most of the broken contacts failed in tension, which was consistent with the theoretical expectations [46]. It is noteworthy that there was a circular area under the loading plane with a significantly smaller number of cracks, which was similar to the result of the experiment (see Figure 6a). Figure 14b and 16c show enlarged images of the cone core after the overall breakage to study the micromechanical process and mechanism underlying the cone core formation. The small images shown in Figure 14 and 16 are zoomed sections of Figure 14a. The transparency of elementary particles was set to 60%, and only particles and cracks in a limited area in front of the cross-section can be observed.
It is clear in Figure 14b and 16c that most of the bonded contacts in the area under the loading plane did not break after the overall breakage. The image acquired from the cross-section perpendicular to the fracture plane also indicates that the shape of the unbroken area was similar to a cone. Another important feature is that shear failures of bonds could be found at the edge of the cone core. This indicates that there was a relative movement trend between the cone core and the rest of the specimen. Figure 15 shows the velocity of the elementary particles in Figure 14c. The arrow is located in the center of the particle pointing to the direction of the velocity, and the length of the arrow is linearly related to the velocity of the particle. It is clear that particles in the blue region basically moved downward as a whole. The result indicates that most of the broken contacts failed in tension, which was consistent with the theoretical expectations [46]. It is noteworthy that there was a circular area under the loading plane with a significantly smaller number of cracks, which was similar to the result of the experiment (see Figure 6a). Figure 14b and Figure 16c show enlarged images of the cone core after the overall breakage to study the micromechanical process and mechanism underlying the cone core formation. The small images shown in Figure 14 and Figure 16 are zoomed sections of Figure 14a. The transparency of elementary particles was set to 60%, and only particles and cracks in a limited area in front of the cross-section can be observed.
It is clear in Figure 14b and Figure 16c that most of the bonded contacts in the area under the loading plane did not break after the overall breakage. The image acquired from the cross-section perpendicular to the fracture plane also indicates that the shape of the unbroken area was similar to a cone. Another important feature is that shear failures of bonds could be found at the edge of the cone core. This indicates that there was a relative movement trend between the cone core and the rest of the specimen. Figure 15 shows the velocity of the elementary particles in Figure 14c. The arrow is located in the center of the particle pointing to the direction of the velocity, and the length of the arrow is linearly related to the velocity of the particle. It is clear that particles in the blue region basically moved downward as a whole. Appl. Sci. 2019, 9, x FOR PEER REVIEW 13 of 16 Figure 15. Force-displacement curve in numerical simulation The bond states at different loading stages are shown in Figure 16. The elementary particles and tensile failures on the cross-section B-B in Figure 14a are displayed. All shear failures are shown. The results show that local breakages first occurred at the surface and inside the specimen under the loading plate; most of the cracks exhibited tensile failure, while some cracks in the annular area centered around the loading point exhibited shear failure (see Figure 16a). Since the load continued to increase, shear cracks developed to the inner of the specimen and eventually formed a cone core (see Figure 16b and 18c). The cone core moved downward under the loading plate and ultimately split the specimen into pieces.

Discussion
Through the experiment, the contact breakage behavior of gypsum particles was analyzed. The force-displacement curve was almost linear, and the breakage was an instantaneous brittle fracture. A significant feature of particle contact breakage was the observation of a cone core formed at the contact point. Limited by the experimental conditions, the formation process and mechanism of the cone core could not be observed by conventional means. The results show that local breakages first occurred at the surface and inside the specimen under the loading plate; most of the cracks exhibited tensile failure, while some cracks in the annular area centered around the loading point exhibited shear failure (see Figure 16a). Since the load continued to increase, shear cracks developed to the inner of the specimen and eventually formed a cone core (see Figure 16b and 18c). The cone core moved downward under the loading plate and ultimately split the specimen into pieces. The bond states at different loading stages are shown in Figure 16. The elementary particles and tensile failures on the cross-section B-B in Figure 14a are displayed. All shear failures are shown. The results show that local breakages first occurred at the surface and inside the specimen under the loading plate; most of the cracks exhibited tensile failure, while some cracks in the annular area centered around the loading point exhibited shear failure (see Figure 16a). Since the load continued to increase, shear cracks developed to the inner of the specimen and eventually formed a cone core (see Figure 16b and 18c). The cone core moved downward under the loading plate and ultimately split the specimen into pieces.

Discussion
Through the experiment, the contact breakage behavior of gypsum particles was analyzed. The force-displacement curve was almost linear, and the breakage was an instantaneous brittle fracture. A significant feature of particle contact breakage was the observation of a cone core formed at the contact point. Limited by the experimental conditions, the formation process and mechanism of the cone core could not be observed by conventional means.

Discussion
Through the experiment, the contact breakage behavior of gypsum particles was analyzed. The force-displacement curve was almost linear, and the breakage was an instantaneous brittle fracture. A significant feature of particle contact breakage was the observation of a cone core formed at the contact point. Limited by the experimental conditions, the formation process and mechanism of the cone core could not be observed by conventional means.
The results of the numerical simulation now provide evidence for the mechanism of the formation of cone cores. The process of contact breakage for gypsum particles under ball-plane contact conditions can be described as follows: Due to the stress concentration, local breakage occurred at the point where the specimen was in contact with the loading plate (see Figure 17b). The breaking part of the specimen was confined by the constraint of the surrounding specimen and was gradually compressed into the specimen during the loading process. The breaking part pressed into the inner part continually compressed the surrounding specimen. The shear cracks developed inside the specimen and gradually formed a conical region (see Figure 17c). The cone core was continuously squeezed downwards into the specimen, and a tensile force perpendicular to the loading direction was exerted on the specimen inside. When the tensile stress inside the specimen exceeded its strength, the fracture formed and rapidly expanded. Finally, overall breakage occurred (see Figure 17d). The results of the numerical simulation now provide evidence for the mechanism of the formation of cone cores. The process of contact breakage for gypsum particles under ball-plane contact conditions can be described as follows: Due to the stress concentration, local breakage occurred at the point where the specimen was in contact with the loading plate (see Figure 17b). The breaking part of the specimen was confined by the constraint of the surrounding specimen and was gradually compressed into the specimen during the loading process. The breaking part pressed into the inner part continually compressed the surrounding specimen. The shear cracks developed inside the specimen and gradually formed a conical region (see Figure 17c). The cone core was continuously squeezed downwards into the specimen, and a tensile force perpendicular to the loading direction was exerted on the specimen inside. When the tensile stress inside the specimen exceeded its strength, the fracture formed and rapidly expanded. Finally, overall breakage occurred (see Figure 17d). During the testing process, the position of the final fracture is affected by the original cracks or defects existing in the specimen. Under the action of tensile stress, a fracture forms at the place with the lowest strength.
Since the total number of elementary particles is limited, the simulation of local breakage is still not sufficiently precise. However, the trade-off is acceptable given the huge amount of computation. In future work, we will generalize this model to other brittle rock materials to further investigate the mechanism of particle breakage.

Conclusions
In this paper, the mechanism of spherical gypsum particle breakage under ball-plane contact condition was experimentally and numerically studied.
Particle contact tests of high-strength gypsum were carried out to investigate the characteristics of particle contact breakage under ball-plane conditions. Unconfined compression tests and splitting tests were performed to determine the physical parameters of the test material. The corresponding numerical simulation was performed using the DEM to study particle contact breakage in a microscopic way. The numerical specimen adopted a flat-joint model to simulate gypsum material. The validity of the methodology was verified by modeling the unconfined compression test and splitting test. The numerical model fitted with the results of the standard rock tests in several aspects, and this feature was also observed during the modeling of the particle contact test.
The conclusions are: (1) The particle breakage a under ball-plane contact condition was an instantaneous brittle fracture. A cone core could be found under the contact area after breakage.
(2) The flat-joint model in PFC is suitable for simulating brittle rock material, such that the results of simulation were well fitted with the results of the experiment. By means of numerical simulation, the mechanism of particle contact breakage was observed from a micromechanics angle. During the testing process, the position of the final fracture is affected by the original cracks or defects existing in the specimen. Under the action of tensile stress, a fracture forms at the place with the lowest strength.
Since the total number of elementary particles is limited, the simulation of local breakage is still not sufficiently precise. However, the trade-off is acceptable given the huge amount of computation. In future work, we will generalize this model to other brittle rock materials to further investigate the mechanism of particle breakage.

Conclusions
In this paper, the mechanism of spherical gypsum particle breakage under ball-plane contact condition was experimentally and numerically studied.
Particle contact tests of high-strength gypsum were carried out to investigate the characteristics of particle contact breakage under ball-plane conditions. Unconfined compression tests and splitting tests were performed to determine the physical parameters of the test material. The corresponding numerical simulation was performed using the DEM to study particle contact breakage in a microscopic way. The numerical specimen adopted a flat-joint model to simulate gypsum material. The validity of the methodology was verified by modeling the unconfined compression test and splitting test. The numerical model fitted with the results of the standard rock tests in several aspects, and this feature was also observed during the modeling of the particle contact test.
The conclusions are: (1) The particle breakage a under ball-plane contact condition was an instantaneous brittle fracture. A cone core could be found under the contact area after breakage. (2) The flat-joint model in PFC is suitable for simulating brittle rock material, such that the results of simulation were well fitted with the results of the experiment. By means of numerical simulation, the mechanism of particle contact breakage was observed from a micromechanics angle. (3) During the loading process, the cone core would be formed by shear cracks under the contact region. The tensile stress caused by the loading force and by the penetration of the cone core would finally exceeds its strength, that the overall breakage occurs.