Introduction

Coarse-grained soils, widely utilized in critical infrastructure projects including earth-rock dams, highways, and bridge piers1,2, demonstrate exceptional engineering properties such as high compaction density, substantial shear strength, and superior bearing capacity. Researchers worldwide have employed experiments1,2,3 and numerical simulations4,5,6,7,8 to rigorously investigate the strength and deformation characteristics of coarse-grained soils since the 1980s. Initially, scholars conducted research on coarse-grained soil using traditional methods. However, they discovered significant differences in the properties of coarse-grained soil compared to ordinary soil. The primary reason for this is that particle breakage has a considerable impact on the mechanical characteristics of coarse-grained soil9,10,11,12,13. After discovering the particle breakage characteristics, most scholars adopted macro-research methods to analyze the impact of particle breakage on the mechanical properties of coarse-grained soil. However, studies conducted solely from a macro perspective have encountered a bottleneck.

In recent years, some scholars have initiated research into particle breakage from a particle-scale perspective in search of a breakthrough14,15,16. The main experimental methods include the triaxial test17,18, the shear test19,20,21, the uniaxial compression test22, the particle contact test23,24,25,26,27,28,29 and the centrifuge modeling test30. Unlike other experimental methods that focus on soil bodies composed of multiple particles (particle groups) as the research object, the particle contact test is a method used to study the crushing characteristics of a single particle26. The influence of particle shape on contact breakage is highly complex. Therefore, in the initial phase of basic research, it is often the case that researchers frequently discuss spherical particles rather than naturally shaped particles7,23,24,25. Unlike traditional fragmentation tests on individual rock samples subjected to a pair of axial rigid point loads, Zhou25 proposed a contact test method between two coarse-grained soil particles to more accurately simulate the fragmentation process of coarse-grained soil particles in contact during actual engineering scenarios. And this method is being adopted by an increasing number of scholars to observe the evolution of spherical particle crushing. Numerous studies indicate that both contact point compaction and global elastic deformation occur during the process of particle breakage, often resulting in the formation of a conical nucleus. However, exploring the internal mechanical characteristics of the particles during the crushing process solely through experimental means is difficult. Theoretical analysis can facilitate an understanding of the distribution of the internal stress state of the particle during the particle contact test. Given the ubiquity of particle breakage processes in different materials, there is a pressing need to propose a contact model to describe the whole process of particle breakage.

The earliest contact model for two spherical particles was the Hertz model31, which considering them as fully elastic and obtained the relationship between pressure and deformation. However, since the particles are fully elastic, discrepancies arise between the predicted stress state at the contact point and the actual conditions. Based on the Hertz model, Hiramatsu and Oka31 provided an analytical solution for an isotropic elastic spherical body compressed diametrically, which was modeled as uniform normal stress acting diametrically on two small areas of the sphere surface. Recognizing that the assumption of using uniform stress for modeling the two applied diametral point loads may lead to either overestimation or underestimation of the maximum tensile stress along the axis of loading, depending on the assumed contact area, some scholars32,33,34,35 optimized the stress state at the contact point based on Hiramatsu and Oka31 model. Particle contact tests demonstrate that the particle breaking process is accompanied by the formation and penetration of the conical nuclei, which is crucial in the mechanism of particle breakage7,28. The aforementioned models improved the boundary conditions to obtain the stress state within the particle near the center of the sphere more accurately, but significant errors persist in solving the stress state in the area near the contact point. Furthermore, the formation process of the conical nucleus includes plastic deformation, which traditional models fail to address effectively. Since the conical nucleus plays a key role in the mechanism of particle breakage, a contact-breakage model considering the formation mechanism of the conical nucleus must be developed.

This paper proposes a contact-breakage (CB) model to delineate the entire process of contact and crushing of spherical coarse-grained soil particles. The CB model is comprised of three phases: the local compaction phase, the elastic deformation phase, and the integral crushing phase, each characterized by distinct mechanical properties. This model takes into account the formation mechanism of the conical nucleus as well as the evolution of boundary conditions during the contact process. Subsequently, the model is validated through test results derived from particle contact tests.

Contact-breakage model

Drawing upon previous research work7,28, the CB model incorporates a three-phase assumption: the local compaction phase, the elastic deformation phase, and the integral crushing phase. During the local compaction phase, the deformation of the particle is primarily attributed to the compaction at the contact point, with overall elastic deformation being negligible. In the elastic deformation phase, the deformation of the entire particle is predominantly elastic, with compaction at the contact point being negligible. During the integral crushing phase, the particle undergoes fragmentation into two pieces, accompanied by the detachment of conical nuclei. A detail description of these three phases is provided in the subsequent section.

Local compaction phase

Numerous factors influence the contact breakage of particles, including particle size, material and shape, rendering it challenging to analyze them all simultaneously. The present study focuses solely on the conditions of contact fracture between two identical spherical particles. Given that both particles are isotropic and possess identical size and material properties, the crushing process of the spheres is symmetrical. Consequently, it suffices to analyze one of the particles. Initially, the place where the two spherical particles come into contact is a point. With the generation and gradual increase in contact pressure, stress concentration arises, leading to the transformation of the contact point into a gradually expanding circular compacted area. As shown in Fig. 1, the overall deformation s of the sphere can be decomposed into the deformation sp resulting from the local compaction and the elastic deformation se of the entire sphere.

Fig. 1
figure 1

A sketch showing the deformation s in the process of contact. Fn is the normal pressure between the two spheres. q is the pressure per unit area. The radius of the circular plane produced by local compaction is a. H is the height of the conical nucleus. R is the radius of the spherical particles.

Full size image

During the phase of local compaction, the authors employ the parameter δ to signify the degree of difficulty associated with local compaction, referring to it as the compaction modulus. The physical significance of this parameter can be articulated through the following equation.

$$W_{c} = \delta V,$$
(1)

where Wc is the input energy for local compaction, δ is the crushing modulus of the material and V is the total locally compacted volume. Due to the presence of pressure at the contact point of two spherical particles, a flat circular contact area forms, expanding as the pressure increases. Consequently, the size of the contact area between the two particles increases, resulting in alterations to the boundary conditions as the loading process progresses.

Initially, the pressure Fn is minimal, rendering the elastic deformation of the sphere negligible. Therefore, the following equation is satisfied.

$$s \approx s_{p} .$$
(2)

During this phase, the work performed by the pressure Fn is converted into the internal energy of the local compaction zone. According to Eq. (1), the following relationship is established:

$$W_{c} = \mathop \int \limits_{0}^{s} F_{n} {\text{d}}s = \delta V,$$
(3)
$$V = \mathop \int \limits_{0}^{s} \pi a^{2} {\text{d}}s = \mathop \int \limits_{0}^{s} \pi \left[ {R^{2} - (R - s)^{2} } \right]{\text{d}}s.$$
(4)

At the very beginning of the local compaction stage, both the contact force Fn and the relative displacement s between the two spheres are zero. As the compaction process progresses, the density and hardness of the compressed region gradually increase, leading to a corresponding gradual increase in the crushing modulus δ of this region, as shown in Fig. 2. By combining Eqs. (3) and (4) and incorporating past research findings to modify the δ value, the following equations are obtained:

$$F_{n} = \frac{\pi }{3}\delta Rs^{3} - \frac{\pi }{12}\delta s^{4} ,$$
(5)
Fig. 2
figure 2

A schematic diagram illustrating the process of partial compaction between two spheres. The shade of color indicates the degree of material compaction; darker colors represent materials that are more densely compressed, corresponding to higher density and hardness.

Full size image

For the purpose of facilitating calculations in this paper, a reasonable assumption is proposed: upon the occurrence of local compaction, a circular planar compaction zone forms at the contact point between the spheres. Within this plane, the normal force exists in the form of a uniformly distributed load, denoted as q, and the uneven distribution of load due to the surface roughness of the spheres is neglected. As a result, the following equation is obtained:

$$q = \frac{{F_{n} }}{{\pi a^{2} }} = \frac{{4\delta Rs^{2} - \delta s^{3} }}{24R - 12s}.$$
(6)

As the normal pressure increases, the compacted area gradually increases. The compaction process at the contact point between two spheres, induced by stress concentration, exhibits a boundary. As this process evolves to a certain extent, it tends to terminate. At any time during this phase, the work required for Fn to compact the area with a thickness of ds is as follows:

$${\text{d}}W_{c} = \pi \delta \left[ {R^{2} - (R - s)^{2} } \right]{\text{d}}s.$$
(7)

It is assumed that the global elastic deformation of ds on a sphere necessitates the work performed by Fn denoted as dWe, i.e., dWe = Fnds. This paper makes the following assumptions: When dWc < dWe, only local compaction occurs on the sphere, with overall elastic deformation being negligible; once dWc ≥ dWe, local compression ceases, and all the work performed by Fn contributes to the overall elastic deformation; additionally, when s = s0, dWc = dWe. Consequently, s0 represents the critical displacement that delineates the transition between the local compaction phase and the elastic deformation phase.

Elastic deformation phase

When s > s0, the phase of elastic deformation initiates, during which the locally compacted area ceases to increases, and all displacement at the contact point arises solely from the elastic deformation of the entire sphere. At this juncture, the force exerted on the contact point can be regarded as a uniform pressure q acting on a circular area with a radius of a. Solutions for the stresses within the half-space along the line normal to and passing through the center of the contact area, as well as the displacement of the contact area’s center, are therefore of particular interest, given that a is significantly smaller than R and the solutions facilitate straightforward mathematical manipulations. Hence, the solution for the central displacement of the compacted area upon entering the elastic deformation phase is as follows:

$$s^{\prime } = \frac{{2\left( {1 - \nu^{2} } \right)}}{E}qa = \frac{{2\left( {1 - \nu^{2} } \right)}}{\pi Ea}F_{n} ,$$
(8)

where E and ν are Young’s modulus and Poisson’s ratio of the material, respectively. Substitute equation \(s = s_{0} + s^{\prime }\) into Eq. (8) and the following equation is obtained:

$$F_{n} = \frac{\pi Ea}{{2\left( {1 - \nu^{2} } \right)}}\left( {s - s_{0} } \right) = \frac{{\pi E\sqrt {R^{2} - \left( {R - s_{0} } \right)^{2} } }}{{2\left( {1 - \nu^{2} } \right)}}\left( {s - s_{0} } \right).$$
(9)

Based on the previous analysis, s0 is the critical displacement separating the local compaction phase and the elastic deformation phase. Consequently, the following equation is obtained:

$$\left. {\frac{{{\text{d}}F_{n} }}{{{\text{d}}s}}} \right|_{{s = s_{0}^{ - } }} = \left. {\frac{{{\text{d}}F_{n} }}{{{\text{d}}s}}} \right|_{{s = s_{0}^{ + } }} .$$
(10)

According to the Eqs. (5) and (9), the following equations are obtained:

$$\frac{{{\text{d}}F_{n} }}{{{\text{d}}s}} = \pi \delta Rs^{2} - \frac{\pi }{3}\delta s^{3} \quad {\text{for}}\;0 < s \le s_{0} \quad {\text{and}}$$
(11)
$$\frac{{{\text{d}}F_{n} }}{{{\text{d}}s}} = \frac{{\pi E\sqrt {R^{2} - \left( {R - s_{0} } \right)^{2} } }}{{2\left( {1 - \nu^{2} } \right)}}\quad {\text{for}}\;s > s_{0} .$$
(12)

Then the following equation is obtained by combining Eqs. (10)–(12):

$$\frac{1}{9}s_{0}^{5} - \frac{2}{3}Rs_{0}^{4} + R^{2} s_{0}^{3} + \frac{{E^{2} s_{0} }}{{\left[ {2\delta \left( {1 - \nu^{2} } \right)} \right]^{2} }} - \frac{{2E^{2} R}}{{\left[ {2\delta \left( {1 - \nu^{2} } \right)} \right]^{2} }} = 0.$$
(13)

Given that s0 is significantly smaller than R, Eq. (13) can be simplified to the following equation:

$$R^{2} s_{0}^{3} + \alpha s_{0} - 2\alpha R = 0,\quad {\text{where}}\;\alpha = \left[ {\frac{E}{{2\delta \left( {1 - \nu^{2} } \right)}}} \right]^{2} .$$
(14)

Equation (14) reveals that α is a dimensionless parameter, that plays a key role in the size of the compacted area (conical nucleus) of particles composed of various materials and is determined by Young’s modulus, Poisson’s ratio and crushing modulus. According to the Girolamo Cardano formula, and since \(\Delta = \left( { - \frac{\alpha }{R}} \right)^{2} + \left( {\frac{\alpha }{{3R^{2} }}} \right)^{3} > 0\), a unique real number solution for s0 is obtained:

$$s_{0} = \sqrt[3]{{\frac{\alpha }{R} + \sqrt {\frac{{\alpha^{2} }}{{R^{2} }} + \frac{{\alpha^{3} }}{{27R^{6} }}} }} + \sqrt[3]{{\frac{\alpha }{R} - \sqrt {\frac{{\alpha^{2} }}{{R^{2} }} + \frac{{\alpha^{3} }}{{27R^{6} }}} }}.$$
(15)

Subsequently, based on the geometric relationship, the area S of the conical nucleus base can be expressed as:

$$S = \pi a^{2} = \pi \left[ {R^{2} - \left( {R - s_{0} } \right)^{2} } \right].$$
(16)

Integral crushing phase

In the elastic deformation phase, the solutions for the stresses within the sphere along the line normal to and through the center of the compacted area are as follows:

$$\begin{aligned} \sigma_{z} & = q\left[ {1 - \frac{{z^{3} }}{{\left( {a^{2} + z^{2} } \right)^{3/2} }}} \right], \\ \sigma_{\rho } & = \sigma_{\varphi } = \frac{q}{2}\left[ {\left( {1 - 2\nu } \right) + \frac{{z^{3} }}{{\left( {a^{2} + z^{2} } \right)^{3/2} }} - \frac{2(1 + \nu )z}{{\left( {a^{2} + z^{2} } \right)^{1/2} }}} \right]. \\ \end{aligned}$$
(17)

σρ, σφ and σρ are normal stress components defined using cylindrical coordinates (see Fig. 3).

Fig. 3
figure 3

Cylindrical coordinates ρ ~ φ ~ z and orthogonal normal stresses σρ, σφ and σz. ρ and z are respectively the horizontal and vertical distance from the center point of the compaction area to the ___location of interest. φ is the horizontal rotation angle of the ___location of interest relative to the x-axis.

Full size image

In this study, we assume the applicability of the simple two-parameter multiaxial failure criterion for brittle materials proposed by Christensen36,37. The criterion stipulates that a material is at failure when the following inequality is satisfied:

$$\frac{\psi K}{{\sqrt 3 }}\sigma_{ii} + \left( {1 + \psi } \right)^{2} s_{ij} s_{ij} > \frac{{K^{2} }}{1 + \psi },$$
(18)

where σii represents the first invariant of the stress tensor and sijsij denotes the second invariant of the deviatoric stress tensor. Furthermore, based on Eq. (17), the following equations are obtained:

$$\begin{aligned} & \sigma_{ii} = 2q\left\{ {\nu - \frac{2\nu + 1}{{2\left[ {1 + (z/a)^{2} } \right]^{1/2} }} - \frac{1}{{2\left[ {1 + (z/a)^{2} } \right]^{3/2} }} + 1} \right\}, \\ & s_{ij} s_{ij} = \frac{{q^{2} }}{3}\left\{ {\nu - \frac{2\nu + 1}{{2\left[ {1 + (z/a)^{2} } \right]^{1/2} }} + \frac{1}{{\left[ {1 + (z/a)^{2} } \right]^{3/2} }} - \frac{1}{2}} \right\}^{2} . \\ \end{aligned}$$
(19)

The parameter ψ and Κ are defined as follows36,37:

$$\begin{aligned} & \psi = \frac{{\left| {\sigma_{c} } \right|}}{{\sigma_{t} }} - 1, \\ & K = \frac{1 + \psi }{{\sqrt 3 }}\left| {\sigma_{c} } \right|, \\ \end{aligned}$$
(20)

where σt and σc represent the uniaxial tensile strength and uniaxial compressive strength of the material, respectively. Russell and Muir34 defined a specific quantity of the available strength characterized by the value of Κ that the multiaxial stress state within the sphere mobilizes as Κmob, and a simple approximate expression for Κmob is as follows:

$$K_{mob} = - \frac{{\sqrt 3 \left( {1 + \psi } \right)^{2} }}{\psi } \frac{{s_{ij} s_{ij} }}{{\sigma_{ii} }}.$$
(21)

In this case, the criterion is that a material is at failure when the following equation is satisfied:

$$K_{mob} = K.$$
(22)

Combining Eqs. (19)–(21), the plot of Κmob within the sphere along the line normal to and through the center of the base of the conical nucleus is presented in Fig. 4. The plot reveals that when Κmob at a point on the central axis of the inner sphere, which passes through the center of the bottom surface of the conical nucleus, reaches Κ, the entire sphere undergoes crushing. This point is where the crystal in the material begins to fracture and can also be considered as the apex of the conical nucleus. Suppose that Κmob is at a maximum when the following equation is satisfied:

$$z = \xi a.$$
(23)
Fig. 4
figure 4

A plot of Κmob along the line normal to and through the center of the loaded area for σc = 50 MPa, σt = 3 MPa, ν = 0.3, a = 2.7 mm and R = 25 mm. The position where Κmob = Κ when crushing occurs is the vertex of the conical nucleus.

Full size image

Simultaneously, the equation below is obtained:

$$\xi = \frac{{\sqrt {\left( {8\nu^{2} + 40\nu + 50} \right)\left( {20\nu + 109 - 3\sqrt {309 + 120\nu } - 8\nu^{2} } \right)} }}{{8\nu^{2} + 40\nu + 50}}.$$
(24)

Consequently, the height of the conical nucleus can be determined as follows:

$$H = \xi \sqrt {R^{2} - \left( {R - s_{0} } \right)^{2} } .$$
(25)

ξ varies between approximately 1.03 and 0.94 for values of ν that vary between 0.1 and 0.4. Thus, it can be reasonably assumed that ξ = 1 due to the weak dependence on ν.

Basic equations

Based on the theoretical analysis of each phase, the relationship between the contact force and particle deformation during the particle contact process is expressed in the CB model by the following equation:

$$F_{n} = \left\{ {\begin{array}{*{20}l} {\frac{\pi }{3}\delta Rs^{3} - \frac{\pi }{12}\delta s^{4} , } \hfill & {\quad \left( {0 < s \le s_{0} } \right)} \hfill \\ {\frac{{\pi E\sqrt {R^{2} - \left( {R - s_{0} } \right)^{2} } }}{{2\left( {1 - \nu^{2} } \right)}}\left( {s - s_{0} } \right).} \hfill & {\quad \left( {s > s_{0} \;{\text{and}}\;F_{n} \le F_{cr} } \right)} \hfill \\ \end{array} } \right.$$
(26)

Fcr is the integral crushing force that is determined by Eq. (22). The bottom area and height of the conical nucleus are derived from Eqs. (16) and (25), respectively.

According to the basic equations of CB model, Fig. 5 illustrates the normal force Fn versus the displacement of the central contact point s throughout the entire breakage process of particles composed of different materials. Three distinct phases are evident: local compaction, elastic deformation and integral crushing. The figure reveals that the local compaction phase is an upper concave curve, whereas the curve becomes straight after entering the phase of elastic deformation. The abrupt decline in the curves signifies the integral crushing stage of the particles.

Fig. 5
figure 5

The predicted curve for Fn and s throughout the entire crushing process for spheres composed of limestone, marble, sandstone and gypsum, with diameter 50 mm. The mechanical parameters of these materials are derived from Zhou38. The circles are the boundary point between the local compaction phase and the elastic deformation phase, while the triangles indicate the occurrence of integral crushing. The dotted lines represent post-peak curves of each material that are impossible to predict theoretically.

Full size image

This section elaborates on three phases of the CB model in detail. This model can not only solve the force–displacement relationship in particle contact process, but also solve the conical nucleus size and the breaking force simultaneously. In the next section, particle contact tests are conducted to validate the model.

Model verification

To verify the applicability of the model, a series of particle contact tests were conducted, encompassing a total of 15 groups of experiments. Initially, 5 groups of tests were performed using high-strength gypsum spheres under ball–ball normal contact conditions. The choice of gypsum was motivated by its simplicity in sample preparation, the elimination of natural cracks inherent in rock materials, and its ability to simulate an isotropic brittle material effectively. Subsequently, to further validate the model’s robustness and applicability across different materials, an additional 5 groups of tests were carried out using limestone spheres and another 5 groups using marble spheres, both under similar ball–ball normal contact conditions. These additions allowed for a comprehensive comparison of the model’s predictions with experimental outcomes across a broader range of material properties. The test results from all 15 groups were then compared with predictions based on the CB model, the Hertz model, and the Hiramatsu–Oka model, respectively. Details of the experiments are provided below.

Particle contact tests of high-strength gypsum particles

Materials and methods

The particular gypsum utilized in this study is Sanxin dental plaster produced by Yuyao Xinshi Gypsum Products Co. Ltd, China. At a water content of 25% (calculated as the mass of water divided by the mass of dental stone powder), the compressive strengths of the plaster, after 1 h and in the long term, are approximately 35 MPa and 75 MPa, respectively. During this experiment, gypsum was prepared with a water content of 35% during mixing, and the mixing duration and curing period were 5 min and 14 days, respectively. Excessive particle size may lead to temperature cracking within the specimen due to the exothermic hardening process of gypsum. Particles of insufficient size may hinder the monitoring and observation of the crushing process. Therefore, spheres with a diameter of approximately 48 mm were prepared with a rubber mold, as illustrated in Fig. 6b. Preparation process of gypsum spheres is as follows: (1) Mix water and gypsum powder at a water-to-solid ratio of 35%, and stir thoroughly for 2 min until no bubbles are generated; (2) Pour the mixed slurry into silicone molds while continuously vibrating to eliminate internal bubbles; (3) After sufficient vibration and no bubbles emerging, seal the pouring port to prevent impurities from entering; (4) Allow the mold to rest for 30 min before demolding; (5) Cure the demolded specimens at room temperature (20 °C) for 14 days; (6) After curing, use 1200-grit sandpaper to polish the seams formed during demolding. The fundamental mechanical parameters obtained through rock testing of the high-strength gypsum material utilized in this study are presented in Table 1. Regarding these mechanical parameters, Young’s modulus, Poisson ratio and unconfined compressive strength were determined by compressing solid circular cylinders of the materials with dimensions of 100 mm in length and 50 mm in diameter. Tensile strength was determined through the Brazilian splitting test, using solid circular cylinders of the same material, with dimensions of 50 mm in both length and diameter. To measure the crushing modulus of the gypsum spheres, particle contact tests under ball–plane contact conditions were performed. The loading was terminated when the pressure reached 800 N, and during this process, the work done by the applied pressure and the final compaction volume were recorded. All results are the means of 5 repeated tests.

Fig. 6
figure 6

Some of the apparatuses used were as follows: (a) Rock rheological testing system; (b) Mold for making spheres; (c) Spheres for the test and the arc-shaped supports for fixing them.

Full size image
Table 1 Physical parameters of gypsum.
Full size table

All contact tests were performed utilizing a rock rheological testing system (see Fig. 6a), which can provide a maximum vertical pressure of 500 kN and a horizontal pressure of 300 kN. The minimum controllable deformation value of this rock rheological testing system is 0.001 mm and the minimum controllable mechanical value is 5 N. The system incorporates a high-speed static force–displacement acquisition system, with a maximum acquisition frequency of 100 times per second. In the present study, the acquisition frequency was configured to 5 times per second. The spherical joint support is specially tailored to prevent eccentric compression during the test.

In the test, the upper and lower spheres were fixed with circular arc-shaped supports (see Fig. 6c), and the upper loading plate was controlled to descend at a velocity of 0.001 mm/s until failure occurred in one of the spheres. Before the tests started, the loading plates and fixture were lubricated with oil to minimize friction effects. During the test, the displacement and force applied by the upper loading plate were meticulously recorded. The size of the conical nucleus was measured using a Vernier caliper after the test. A comprehensive discussion of the test results is presented in the subsequent section.

Results and comparison

Prior to conducting the formal test, the experimental results were anticipated based on the CB model. According to Eq. (14), s0 = 0.18 mm, indicating that the loading process transitions into the global deformation phase when the displacement surpasses 0.18 mm. Subsequently, the dimensions of the conical nucleus were derived using Eqs. (16), (24)–(25), yielding S = 26.91 mm2 and H = 2.89 mm. Furthermore, by integrating equations (19–22), we determined that the normal force at crushing is 8.021 kN. Consequently, in Fig. 7, the predicted force–displacement curve (depicted in green) is plotted in accordance with Eqs. (5) and (9). Additionally, the predicted force–displacement curves based on the Hertz model (depicted in blue) and the Hiramatsu–Oka model (depicted in violet) are also plotted. The orange curve represents one of the 5 groups of experimental data obtained from the test, which aligns well with the predicted curve of the CB model. Consistent with the result of the theoretical analysis, the force–displacement curve measured through experimentation demonstrates that the particle breakage process encompasses three phases. In the initial contact process, the relationship between the normal force Fn and displacement s is represented by an upward concave curve. During this stage, local compaction arises at the contact point. Upon entering the elastic deformation phase, local compaction becomes less evident, while global elastic deformation of the sphere becomes dominant, and the curve transitions into an inclined straight line. As the normal contact force increases, the entire sphere undergoes crushing, resulting in a steep decline in the curve. It is evident from Fig. 7 that the CB model exhibits greater consistency with the particle contact test compared to the traditional elastic theory models.

Fig. 7
figure 7

Diagram of normal force Fn and displacement s: predicted curve based on the CB model, Hertz model and Hiramatsu–Oka model and a set of measured data from the test. The circles are the boundary point between the local compaction phase and the elastic deformation phase, while the triangles indicate the occurrence of integral crushing. The dotted lines represent post-peak curves that are impossible to predict theoretically.

Full size image

Particle contact tests of rock particles

To further validate the model’s robustness and applicability across different materials, an additional 10 sets of particle contact tests were conducted, with 5 sets using limestone spherical particles and the other 5 sets using marble spherical particles, both under similar ball–ball normal contact conditions. Unlike gypsum spheres, limestone spheres and marble spheres are both processed by stone factories to ensure uniformity in both quality and size of the specimens. The processing procedure for limestone and marble spheres is as follows: (1) Raw blocks were first cut into cubic samples using a diamond saw; (2) The cubes were then mounted on a CNC spherical grinding machine and shaped into spheres through progressive grinding with 80–3000 grit diamond wheels; (3) Final polishing was performed using a 5000 grit diamond paste to achieve surface roughness Ra < 2 μm; (4) Each sphere underwent dimensional verification using a coordinate measuring machine to ensure diameter uniformity within ± 0.05 mm. The entire process was conducted under controlled temperature (20 ± 2 °C) and humidity (50 ± 5% RH) conditions to minimize thermal and moisture-induced dimensional changes. The limestone used in this experiment exhibits an unconfined compressive strength of 89.6 MPa, a Young’s modulus of 55.3 GPa, a Poisson’s ratio of 0.28, a tensile strength of 5.4 MPa, and a crushing modulus of 43.7 GPa. The corresponding parameters for marble are 70.4 MPa, 39.2 GPa, 0.24, 4.5 MPa, and 35.2 GPa, respectively. The experimental equipment and procedure are identical to those used in the contact test with gypsum sphere particles. Figure 8 compares the relationship between normal force Fn and displacement s predicted by several models with experimentally measured data. It can be observed that for limestone and marble materials, the CB model demonstrates better predictive capability for particle contact tests compared to traditional elastic theory models.

Fig. 8
figure 8

Diagram of normal force Fn and displacement s: predicted curve based on the CB model, Hertz model and Hiramatsu–Oka model and a set of measured data from the test: (a) Particle contact test on limestone spheres; (b) Particle contact test on marble spheres. The circles are the boundary point between the local compaction phase and the elastic deformation phase, while the triangles indicate the occurrence of integral crushing. The dotted lines represent post-peak curves of each material that are impossible to predict theoretically.

Full size image

Table 2 shows the comparison between the measured data and predicted values by the theoretical analysis. For the bottom area and height of the conical nucleus as well as the normal force at crushing, the mean measured and predicted values exhibit relative errors of 9.0%, 6.9%, and 3.0% respectively for high-strength gypsum. For limestone, the corresponding relative errors are 6.8%, 8%, and 7.7%, respectively. Similarly, for marble, the relative errors are 5.6%, 7.7%, and 6.9%, respectively. These results indicate that the theoretical analysis of this paper has good reliability for the calculation of contact breakage between spheres.

Table 2 The means of the bottom area and height of the conical nucleus as well as the normal force at the time of crushing of the tests are listed, with their predicted values listed for comparison.
Full size table

In this section, the authors carried out particle contact tests using high-strength gypsum spheres, limestone spheres and marble spheres to simulate an isotropic material with a continuous medium. The experimental data were compared with the theoretical predictions in terms of the force–displacement curve, size of the conical nucleus and the crushing force. The experimental results matched the theoretical predictions very well.

Discussion

In this paper, a CB model for spherical coarse-grained soil particles is presented and validated through particle contact test conducted in “Model verification” section. This model provides predictions for the conical nucleus, the contact force–deformation curve and the failure force. In this section, the basic assumptions, applicability and limitations are discussed.

Classical contact models based on the Hertz model assume that particles are completely elastic and, therefore, fail to describe the compaction process at the contact point, as local compaction involves plastic deformation. Consequently, these models yield significant errors in calculating the stress state near the contact point. To address the aforementioned problem, the CB model assumes that the external input energy required to compact a unit volume remains constant as the contact force increases, as formulated in Eq. (27). The reason is that the area of the compaction zone at the contact point is substantially smaller than the particle size.

$$\delta = \frac{{{\text{d}}W}}{{{\text{d}}V}} = \frac{W}{V}.$$
(27)

In the CB model, the three phases are categorized based on the deformation type that dominates in each phase. Indeed, both plastic deformation at the contact point and elastic deformation of the entire particle coexist before breakage. During the local compaction phase, the elastic deformation of the particle is substantially smaller than the plastic deformation at the contact point. To facilitate calculation, the model disregards the elastic deformation of the particle in this phase. According to Eq. (7), the energy required for compaction per unit thickness increases markedly as particle contact develops. Furthermore, the crushing modulus increases significantly when the local compaction process reaches a certain stage, further enhancing the difficulty of local compaction. Consequently, in the elastic deformation phase, the CB model neglects the compaction at the contact point. When particle crushing occurs, cracks undergo three stages: generation, development and penetration, as illustrated in Fig. 9. This paper focuses on brittle materials as the research object, with a quasi-static loading rate. Consequently, the time from crack generation to penetration is substantially shorter than the loading time. The CB model disregards the development process of cracks and asserts that the integral crushing occurs once Eq. (22) is fulfilled.

Fig. 9
figure 9

The development process of cracks when breakage occurs.

Full size image

In the local compaction phase, the contact mode of two particles changes from point contact to surface contact (see Fig. 10). To facilitate the determination of the internal stress state of the particle, the CB model assumes a uniform distribution of the load across the contact surface of the two particles.

Fig. 10
figure 10

The variation of boundary conditions.

Full size image

The CB model is subject to certain conditions for use. It is exclusively applicable to spherical particles and requires further refinement for particles of other shapes. However, the piecewise solution approach embodied in the model can be adapted to the particle breakage of any shape under point contact conditions (see Fig. 11). The model can be effectively applied to two spherical particles of identical size, while further enhancements are necessary for particles of varying sizes. In addition, the model’s applicability is limited to brittle and homogeneous rock materials, resulting in significant errors for non-brittle rocks or rocks with numerous natural fractures.

Fig. 11
figure 11

A diagram of particle breakage of various shapes under point contact conditions.

Full size image

Conclusion

This study presents a comprehensive contact-breakage (CB) model for two brittle spherical coarse-grained soil particles, considering the formation mechanism of the conical nucleus. The main conclusions are summarized as follows:

  1. 1.

    The proposed CB model successfully delineates the entire process of spherical particle breakage into three distinct phases (local compaction, elastic deformation, and integral crushing), incorporating the formation mechanism of the conical nucleus and the evolution of boundary conditions during contact, which provides a more comprehensive understanding of particle breakage compared to traditional elastic theory models.

  2. 2.

    Experimental validation using high-strength gypsum, limestone, and marble spheres demonstrates the model’s high accuracy in predicting the size of the conical nucleus, the contact force–deformation curve, and the force at failure, with the CB model outperforming traditional elastic theory models particularly in describing the stress state near the contact point.

  3. 3.

    The integration of the conical nucleus formation mechanism into the particle contact-breakage model provides a new perspective for understanding and predicting particle breakage in coarse-grained soil particles, while the model’s current applicability is limited to spherical particles of identical size and brittle, homogeneous rock materials.

  4. 4.

    This study contributes to advancing particle breakage research by enhancing the understanding of its mechanisms and offering potential insights for engineering applications with coarse-grained soils, with future research directions including extending the model to particles of varying shapes and sizes, non-brittle rocks, and rocks with natural fractures.