## Abstract

The particle size plays an important role in the mechanical behavior of granular geomaterials. With respect to the crushable granular soils, the particle breakage under different pressures leads to the variation of the particle size distribution of soils. Therefore, it is essential to accurately model the particle fracture process to reproduce the key response of granular soils in the simulation. In this study, two simulated approaches, cluster method and replacement method, were performed to investigate the particle size effect for the crushable granular soils by three-dimensional (3D) discrete element method (DEM). DEM simulations were validated through the available results of diametric and oedometric compression tests with the help of the Weibull statistics. Based on the evaluation at particle and assemble-scale, it is found that both the cluster method and replacement method could simulate the particle strength and size dependence of sample failure with a good agreement with experimental results. The particle breakage phenomenon is also successfully modeled in DEM simulation and analyzed by the breakage index. Meanwhile, the computational efficiency was analyzed compared with two simulated methods. The findings in this study can provide a viable approach to investigate the particle breakage effect of the crushable geomaterials.

## 1. Introduction

It has been shown that the mechanical behavior of granular materials is closely related to the characteristic of individual particles, including particle size, particle shape, particle breakage, and particle arrangement [1]. The breakage of particles and the reduction of the particle size are the main mechanism dominating the constitutive behavior of crushable granular assemblies [2, 3]. These microscale evolutions are of significant importance for their engineering applications, such as railway tracks [4], rockfill dams [5], embankments [6], displacement piles (Liu et al., [7, 8]), landslides ([9]; Liu et al., [7, 8]), and oil/gas recovery [10]. The particle size of the granular material plays a crucial role in determining the macroscale response of polydisperse assemblies, which is called the particle size effect.

A lot of studies have been conducted on the particle size effect of granular materials over the past decades because it is essential to consider this phenomenon to provide useful engineering solutions. Most of the experimental tests were carried out through the diametric and oedometric compression experiments. Nakata et al. [11] performed a series of diametric compression tests on quartz and feldspar samples with particle diameter ranging from 0.85 mm to 2.0 mm to evaluate the response of a single particle during crushing. The results demonstrated that (1) the Weibull prediction shows good agreement with the quartz data of the variation of strength with particle size and (2) the characteristic crushing strength $\sigma nf$ decreases with the particle size $d$. Similar results were also reported by Lee [12], who carried out a series of diametric compression tests to characterize the relationship between mean tensile strength and particle size on the limestone and sand. McDowell and Amon [13] presented an analysis of the Weibull statistics applied to tensile failure of soil particles compressed between flat platens as well. The results indicated that 37% of compression strength $\sigma co$ decreases with the increasing particle size. In terms of the particle assemblies, Nakata et al. [14] conducted one-dimensional (1D) compression tests on silica sand with median particle diameter $d50$ ranging from 0.3 mm to 1.7 mm, and they found that the yield stress $\sigma vy$ and the stress at maximum compression index $\sigma vcc=max$ reduce with the increase of median particle size $d50$. The same results were concluded by McDowell [15], who summarized the yield stress for each aggregate with different initial gradings from 0.3 mm to 2 mm in 1D compression test.

The experimental test is regarded as a reasonable approach in studying mechanical response of particle fracture. However, it is hard to replicate the particle shape, packing, and stress condition, thus deteriorating the predictive capability of models [16]. The difficulty in preparing homogeneous samples can be averted by using the discrete element method (DEM), which has been one of the most common numerical approaches to investigate the macro- to microbehavior of granular materials [17]. Several ways to simulate the particle breakage have been applied in numerical investigation up to now, which includes (1) the crushable particles represented by a cluster of bonded element spheres ([18–21]), (2) the breakage particles replaced by predefined smaller fragments based on the preset failure criterion ([22–25]), (3) the combined FDEM modeling method with nonthickness cohesive interface elements inserted between the finite element of each particle ([26]), and (4) the failure of bonds modeled by a damage model using peridynamics [27]. The particle breakage has been extensively studied in recent years by using DEM in the powder technology area. The process of grain crushing in granular media is featured by accompanied evolutions of both grain size and shape. Take natural sands as an example, it has been well accepted that continuous particle crushing would lead to particle size evolving towards a fractal, in other words self-similar, distribution [28]. The grain-scale strength is dependent on particle size as a result of statistical distribution of microdefects in the bulk of grain. Zhu and Zhao [28] investigated the particle shape evolution and its interplay with the crushing characteristics during continuous crushing of the granular sand by employing a hybrid peridynamics and nonsmooth contact dynamics approach to simulate particle breakage. Shi et al. [29] investigated the effects of intermediate principal stress ratio $b$ and particle breakage on the critical state behaviors of granular materials using the octahedral shear stress crushing criterion. However, few attention has been drawn to the particle size effect, by employing these breakage models.

For particle breakage modeling with DEM in previous studies, few attention has been drawn to particle size effect [30]. Wang and Arson [2] presented that the particle size effect of crushable sand particles can be successfully captured by the cluster method with a specific value of ratio, which is defined as element $spheres\u2009size/cluster\u2009size$, smaller than 0.014. Zhou et al. [31] conducted a series of triaxial tests to investigate the size effect of granular materials with the scale of individual particles and assemblies scale using DEM. They found that the particle size effect in peak deviatoric stress and volumetric response does not appear under the low confining pressure, while the higher peak deviatoric stress and lower volume contraction appears in the sample with a smaller particle under high confining pressure. Cil et al. [16] presented a novel modeling strategy of cluster method to study the particle size effect of granular materials in diametric and oedometric compression tests. This alternative strategy relied on simulation of the breakage process at the scale of single particles and assemblies. To sum up, although different materials and test methods were used, they share a common conclusion: the crushing strength decreases with the increasing particle size in diametric compression test; besides, the yield stress of coarse assembly is smaller than that of fine assembly in oedometric compression test. However, to the authors’ knowledge, few studies have been conducted to evaluate particle size effect on the mechanical behavior of granular soils through micromechanical investigation.

The current approaches on simulating brittle process through DEM can be characterized as follows: (1) cluster method, the number of elementary spheres increases with agglomerate size, leading to a low computational efficiency; (2) the replacement method, mass conservation is unable to be guaranteed after the particle fracture, although the computational efficiency is higher. Two important challenges appear when using spheres in association with a particle replacement scheme in DEM [32]. The first is associated with the lack of volume conservation, and the second is the massive superposition of the fragments since spherical fragments cannot properly fit inside the mother particle. Several recent works deal with the mass loss problem. Tavares and Anderson [33] presented a stochastic particle replacement strategy for simulating breakage in DEM. A “dummy” size progeny particle is defined to mimic the unresolved fines from the simulation, thus preventing breakage from creating a mass loss. This breakage model was demonstrated in simulating the breakage of particles in a drop weight tester, whose limitation is that the individual particle families are likely to have a larger volume gain or loss when using a relatively small value of $m$ (i.e., the desired number of parent particles or families). As a summary, despite several novel approaches being proposed to deal with the mass loss during particle breakage simulation, the challenges of keeping mass conservation during replacement still exist. Up to now, no consensus has been made upon the optimum method for modeling particle fracture through DEM.

In this paper, the particle size effect in mechanical behavior of granular soils in diametric and oedometric tests is investigated using the DEM simulation at both macroscopic and microscopic scales. The comparative study is conducted on both cluster and replacement methods for modeling particle crushing from particles to assembly scale. To achieve this goal, the force-displacement curves obtained in diametric compression tests of both methods and experiments were used to calibrate DEM parameters. Based on the calibrated DEM parameters for both methods, the variable stress response in diametric simulation of both methods were plotted against the Weibull statistics. The size-dependent response in oedometric tests was captured by both methods and compared with experimental data.

## 2. DEM Simulations

### 2.1. Brief Introduction of DEM and Contact Model

In this study, the DEM is used to analyze the particle size effect of granular soils in diametric and oedometric tests. The DEM was initially proposed by Cundall and Strack in the 1970s to simulate rock masses and now has been widely used to simulate the behavior of powders, soils, and rocks. All numerical simulations were performed through the commercial DEM software PFC version 6.0 [34], which simulates the mechanical behavior of unbonded particles based on the DEM technique introduced by Cundall and Strack [17]. In the DEM simulations, the interaction between particles is determined by contact, and movement is governed by Newton’s second law. The rigid spheres and rigid walls are allowed to overlap at contacts during the numerical simulations. Contact surface areas are small compared to the surface area of spherical elements. Because the spheres are assumed to be rigid and unbreakable, the fracture of sand particles is here modeled by treating the sand grains as crushable agglomerates composed of bonded elementary spheres by the cluster method [35], termed as Method 1. The replacement method is also used here to simulate the breakage process of sand particles by replacing the breakage grains into a predefined smaller cluster based on the preset failure criterion [23], termed as Method 2. Three constitutive laws determine the contact behavior: a stiffness law, a slip law, and a bonding law. The stiffness model provides a relationship between a contact force and a relative displacement. The slip model characterizes the maximum shear force to support a contact before a slip movement occurs. The bonding model defines the mechanism of bonds, which can undertake force and bending moments. In the simulation, crack propagation is often represented by bond breakage, in which case, bonds disappear after reaching the strength threshold.

The Hertz-Mindlin stiffness model was used in the unbonded group of particles, in which the normal and shear stiffnesses increase with the overlap of spheres. The Hertz-Mindlin slip model was also adopted in unbonded particles. The maximum shear force at contact is defined as $Fmaxs=\mu Fin$, where $\mu $ is the interparticle friction coefficient and $Fin$ is the normal contact force. Therefore, the governing parameters for the unbonded group of particles include the shear modulus $G$, Poisson’s ratio $v$, and the interparticle friction coefficient $\mu $. The parallel bond model was used in the group of bonded spheres, in order to simulate the fracture of sand grains. The bonds among spheres are represented as short beams with a circular cross-section. The bonds can transfer force and bending moment between the interaction particles. The detailed constitutive laws can be referred to Itasca Manual [34] for further details.

### 2.2. Diametric Compression

The validation of the DEM model for diametric tests simulations can be demonstrated qualitatively by comparing the numerical results with the experimental data of single-particle breakage tests by Nakata et al. [11]. The experimental tests were carried out by placing a sand particle between platens and then moving the platen at a constant rate to make the particle crushable. A series of single-particle breakage tests were conducted on quartz according to particle size of 2.0 mm. The compression was applied until the sand grain exhibited catastrophic failure (i.e., the appearance of first fragmentation) followed by a sudden drop in the axial force.

The DEM simulations for the diametric compression tests were carried out here by using two methods to calibrate the parameters in the subsequent modeling. Method 1, termed as the cluster method, is adopted by generating elementary spheres with hexagonal close packing to achieve a very dense packing within a spherical wall container. The spherical wall container has the final size of the agglomerate. The diameter of an elementary sphere is normally chosen in the range of $2Dagg/7$~$Dagg/10$, in which $Dagg$ is the diameter of agglomerate size [39]. In this study, it is taken as $Dagg/7$, whose value can ensure enough elementary ball numbers for each agglomerate and reduce the computation cost. After the generation of elementary balls, the frictionless hexagonal packing granular system is then subjected to consolidation with the gravity of 9.81 m/s^{2} until a state of static equilibrium is achieved. Then, parallel bonds are installed at all existing contacts between the elementary balls whose contact gaps are less than the specified installation gap value $g\xafi$. The value of $g\xafi$ is selected to make sure all the elementary spheres are bonded together. Finally, the spherical container is removed, and the granular soil is subjected to compressive loading between two stiff platens with a constant moving rate. The moving rate is chosen carefully to ensure the soil at a quasistatic status. With respect to Method 2 (i.e., the replacement method), the preparation process of building a DEM model is similar but straightforward compared with that in Method 1. A single spherical ball is compressed between the platen with the same displacement rate to achieve a quasistatic condition. The predefined breakage criteria are checked at a specified cycle interval; once the particle is recognized as breakage, the mother particle will be replaced by several smaller fragments. This breakage model is modified from Wang et al. [40], and the details will be described in Section 2.4.

### 2.3. Oedometric Compression

The stress-strain response of a strain-controlled oedometric compression experiment by Cil et al. [16] is adopted here to validate the behavior of DEM simulations of oedometric compression by both Method 1 and Method 2. The experiments have cylindrical specimens with a diameter of 30 mm and height of 19 mm. Each specimen is subjected to a constant strain rate of 0.1%/min up to axial stress of 100 MPa.

Numerical counterparts of high-pressure oedometric compression tests on three different batches of each sand were performed using the calibrated parameters to examine the link between individual particle failure and macroscopic-yielding characteristics. Initially, an assembly of rigid spherical particles with size distributions similar to those of Q-ROK was randomly generated within a cylindrical vessel (i.e., diameter is $10\xd7Dave$, and height is $5.0\xd7Dave$). Then, after several computational timesteps being imposed to eliminate particle overlaps, in Method 1, each rigid spherical particle was replaced with arbitrarily rotated and scaled crushable aggregates, which are composed of bonded spherical elementary particles. In Method 2, each rigid spherical particle was assigned to the predefined breakage criteria; when the maximum contact force reaches its threshold, the mother particle was replaced by two clusters. Each cluster contains 17 spherical particles, in order to minimize the volume loss during replacement.

The calibrated parameters were assigned after the soil reaches its equilibrium. The quasistatic conditions were granted by imposing a constant strain rate to the top and bottom axial walls. The oedometric compression simulations and experiments are summarized in Table 1. The DEM simulations involve three batches of sand, with varying particle size ranges: Q-ROK #1 (average particle size $Dave$ is 0.35 mm), Q-ROK #2 (average particle size $Dave$ is 0.67 mm), and Q-ROK #4 (average particle size $Dave$ is 1.00 mm). It should be noted that the diameter and height of specimens are set as 10 and 5 times the $D50$ in the simulation. These conditions would have prior evidence to minimize the wall effects. The particle size distribution (PSD) of all sands is plotted in Figure 1.

### 2.4. Parameter Calibration for Breakage Model in DEM

In Method 1, crushable particles are represented by bonded elements, and particle breakage is simulated as the bond breakage. In Method 2, once a particle breaks, it will be replaced by several smaller fragments. The fragment size, shape, and arrangement must be specified in the replacements, which are related to the contact forces and contact positions of the mother particle. Both of the two methods have their advantages: (a) the cluster method can generate fragments of different sizes and shapes but its computational efficiency is low; (b) the replacement method has a higher computational efficiency but requires a reasonable breakage model with suitable criterion.

The breakage model used in this study is presented in Figure 2. This model is proposed from the sequential breakage mechanism characterized from XCT images of oedometer tests and is a combination of cluster method and replacement to simulate particle breakage of multiple generations. In the first generation, once a particle is broken, it is replaced by two clusters with 17 balls bonded with parallel bonds in each cluster. The breakage plane is determined by the directions of principal stresses and the position of the contact with maximum normal force. The 16 smaller spheres are tangent to the sphere in each cluster that represents a fragment which is chosen so as to minimize the mass loss and ensure computational efficiency. Less mass is lost when replacing a particle by its fragments when the cluster contains a large number of small spherical elements. With the cluster of 17 particles, the mass loss is 46%. According to Ciantia et al. [22], it is below the critical mass loss (47%). Sensitivity analysis shows that it is enough to keep mass loss below 47% within the mechanical model to correctly reproduce experimental macroscopic behavior. In the following generations, breakage can be simulated by either replacement method or cluster method depending on the maximum normal contact force and stresses within parallel bonds. As proved by many researchers, the strength of a single particle increases when the particle size decreases [11, 13, 14]. To account for this phenomenon, the Weibull theory is adopted in the breakage model to predict the strength for particles and parallel bonds of various sizes. A detailed description on the failure mechanism and breakage model can be found in McDowell and Amon [13].

The main micromechanical parameters associated with the breakage model are the reference tensile strength in the Weibull theory ($\sigma t0$), the Weibull modulus ($m$), the normal and shear bond strengths ($\sigma \xafc$ and $\tau \xafc$), and the parallel bond normal and shear stiffnesses ($k\xafn$ and $k\xafs$). In this study, $\sigma t0$ represents the tensile strength of the particle with $d=1\u2009mm$, and it ranges from 1.5 to 2.5 MPa to simulate particles with various crushability. Weibull modulus is chosen as 4.5 (Lobo-Guerrero and Vallejo [41]). Both $\sigma \xafc$ and $\tau \xafc$ are assumed to follow the Weibull theory with the same $\sigma t0$ and $m$ with soil particles. The equivalent normal and shear stiffnesses, $\pi R2k\xafn$ and $\pi R2k\xafs$, also equal to the stiffnesses of soil particles. The parameter used in DEM simulation is carefully calibrated, listed in Table 2.

## 3. Results Analysis

### 3.1. Failure of a Particle under Diametric Compression

To assess the potential anisotropy induced by the HCP orientation on mechanical properties of clusters, we compared the results obtained with the set of DEM bond parameters reported in Table 2 for the following orientations: 0°, 30°, 60°, and 90°. The setup of tests for HCP orientation effect is shown in Figure 3. In order to validate the DEM model, the simulated relationships of breakage force and displacement are compared with the measured data in Nakata et al. [11], as shown in Figure 4. In Method 1, different simulated models were considered to investigate the effect of hexagonal close packing (HCP) lattice degrees. It is found that the DEM model with HCP lattice 0° has a close peak breakage force whereas the increase rate of the force is higher than that of the experiment. In comparison, a simulated model with HCP lattice 60° behaves at a close average value of force increase rate with a higher peak breakage force. In this study, the HCP lattice 60° is selected in Method 1 to reproduce the behavior of a particle well. In Method 2, the replaced fragments contain 17 elementary balls and can continue to break until the limitation of minimum ball size. The simulated result of method 2 is similar to that of a simulated model with HCP lattice 0°.

The typical simulation results are shown in Figure 5, in which we can see that due to the symmetry imposed on the model, the particle breaks along the loading axis in both Method 1 and Method 2. When the first peak force is reached, the cluster breaks into four pieces with several smaller fragments by Method 1 and two pieces by Method 2.

The relationships between the normalized characteristic stress and the survival probability ($Ps$) in two methods (Method 1 and Method 2) as well as the experimental data from Nakata et al. [11] are shown in Figure 6. It is obviously found that the survival probabilities ($Ps$) for different agglomerates are almost identical, and there is a good agreement between the experimental data and simulated results when the Weibull statistics is used.

The values of Weibull modulus $m$ can be calculated from the relationship of $lnln1/Ps$ against $ln\sigma $. The simulated results of different agglomerate groups (0.35, 0.67, and 1.0 mm) in Method 1 and Method 2 are plotted in double logarithmic graphs, as shown in Figure 7. It is seen that the $m$ value for crushing stress decreases with particle diameter, and $m$ values obtained from Method 1 are extremely close to those from Method 2 with the exception of 0.67 mm case. From Figure 7, it is also observed that the particle strength has a close relationship with the particle size: the smaller the size, the higher the particle strength.

Based on the proposed function of Equation (3), Figure 8 plots the relationship of characteristic strength and particle size measured in DEM simulations. The slope of the best fit line produced by Method 1 (1.963) is quite close to that obtained by Method 2 (1.994), as shown in Figure 8, which indicates that Equation (3) is also successful in describing the crushable particles.

### 3.2. Macroscopic Yielding

The relationship of stress-strain is seriously concerned for both scholars and engineers. Figure 9 demonstrates the axial stress-strain curves for the DEM simulation (Method 1 and Method 2) and experimental result. The nonlinear behavior of the curve is observed due to the particle crushing. Figure 10 shows the view of samples before and after the axial stress. In general, the simulated results follow the same pattern as the experimental data. For both simulation and experiment, the axial stress-strain curves show a primary linear response, followed by a marked nonlinear change corresponding to the yielding. After yielding, the axial stress-strain curves in the DEM simulations (both Method 1 and Method 2) exhibit obvious fluctuations, which can be attributed to the particle breakage and adjustment. Moreover, it is observed that the yielding region in both the experimental and DEM simulations exhibits variance with the increase of particle size, in which Q-ROK-4 exhibits an earlier yielding response (at around 3 MPa) followed by Q-ROK-2 and Q-ROK-1.

### 3.3. Particle Breakage

The particle breakage would change the original particle size distribution (PSD) of granular soils and results in the variation of internal structure and material properties (compressibility, strength, etc.). Therefore, the particle breakage is closely related to the mechanical behavior of granular soils.

Figure 11 compares the PSDs of experimental data and DEM simulations (cluster method and replacement method) at stresses of 0, 6, 20, 50, and 100 MPa. Generally, there is an acceptable result between simulation and experimental data. It is also shown that particle crushing increases with the applied stress, leading to the percent increase of the fineness of the particle size distribution. When the stress is increased to 100 MPa, in Q-ROK#1 specimen, about 47% of the particles have been broken in the cluster method (Method 1) compared to about 58% and 61% for Q-ROK#2 and Q-ROK#4 specimens, respectively. Similarly, in the replacement method (Method 2), about 45% of the particles have broken for Q-ROK#1 compared to about 52% and 70% for Q-ROK#2 and Q-ROK#4 specimens, respectively. The size distribution in the coarser part varies little in the simulations in comparison with the experiments. This is due to the particle survival: the particle size and the coordination number. Smaller particles are stronger but have fewer contacts, whereas large grains tend to have difficulty to break due to the dominating hydrostatic effects created by their large number of contacts [13, 31, 42]. Large grains have high coordination numbers, which produce well-distributed stresses. The dominating effect of the coordination number over the grain size in determining the probability of grain fracture makes the small grains continually to crush, surrounding the large grains to increase the coordination number [13, 36, 42, 43]. Generally, the PSDs from DEM simulation of Method 2 are consistent with those from DEM simulation of Method 1, which ensures the successful reproduction of the particle breakage phenomena in the compression experiment.

Figure 13 shows the linear relationship of the particle size and the yielding stress of the DEM simulation using both Method 1 and Method 2. The slope values are 1.5 and 1.4 for Method 1 and Method 2, respectively. However, that of experiment with a value of 1.1 is slightly lower than those from DEM simulation. Besides, it shows that particle size has a significant impact on yielding stress: the yielding stress decreases with the particle size, indicating that smaller particles have higher yielding stress.

## 4. Discussions

Basically, there are two types of elastic contact models in DEM: the linear elastic model and the nonlinear elastic model. The linear elastic model can be enlivened as two linear springs providing linear elastic normal (no tension) and shear forces with constant stiffness $kn$ and $ks$. The nonlinear normal elastic model usually refers to the Hertz contact model, which provides a nonlinear contact stiffness. The main difference between the two elastic contact models lies in the normal and shear stiffness. We used the Hertz-Mindlin stiffness model in which the normal and shear stiffness increase with the sphere overlap [34]. The two governing parameters are the shear modulus $G$ and Poisson’s ratio $V$. In the slip model, the maximum shear force at a contact is expressed as $Fmaxs=\mu Fin$, where $u$ is the friction coefficient and $Fin$ is the normal contact force. A detailed description of this model can be referred to in Itasca Consulting Group, Inc. [34].

Particles below a specific size threshold are often excluded from force chains [46] and have a very slight effect on the shear strength when they account for less than 20% total weight [47]. In the present study, the two-generation constraint is set to avoid high computational cost in Method 2, the replacement method. Theoretically, the isolated spherical grains can always break when the breakage criterion is met. However, the large gap between the relative mass of the particles will greatly increase the computational time because the critical time step in DEM is calculated as $tcrit=mmin/k$. With this constraint, the minimum particle diameter in Method 2 is as small as the diameter of the elementary ball in Method 1. These small particles are generally excluded from the force chains ([46]); their breakage is expected to have a very slight effect on the overall behavior of the granular material. However, high computational cost cannot be evaded in this research, especially in Method 1 (i.e., the cluster method). Each agglomerate in the cluster method needs to be filled up with a sufficient number of elementary balls, in order to guarantee the effectiveness and repeatability of the simulation test. As the diameter or number of agglomerates increases in a single test, the number of elementary balls will rise rapidly, resulting in the inevitably high computational cost.

## 5. Conclusions

In this study, we provide a comparative study on the cluster and replacement methods for modeling the size-dependent mechanical behavior of granular soils. A series of diametric compression tests were conducted by controlling the loading velocity of two rigid walls, placed at the top and bottom of the particle. The DEM model of both methods was successfully calibrated through experimental datasets relative to a series of single-particle crushing tests. Based on the comparison and result analysis, main conclusions can be drawn as follows:

- (1)
In diametric compression simulations, the statistical analysis of DEM breakage model reveals that the variability of fracture stress at the beginning of crushing and their size dependence can be successfully characterized by the Weibull statistics

- (2)
With regard to the oedometric compression test, the oedometric response of the two methods was comparable with that of the experiments. The simulation of the mechanical response of particle assemblies shows that both methods captured the size effect in the yielding threshold. The larger particles assemblies exhibit more crushability and have a lower macroscopic-yielding threshold

- (3)
The overall trend of PSD evolution and breakage index of both methods during compression test shows a good agreement with the experiments

- (4)
The computational cost of Method 1 is higher than that of Method 2 in both diametric compression and oedometric compression.

## Data Availability

All data, models, or code that support the findings of this study are available from the corresponding author upon reasonable request.

## Conflicts of Interest

The authors declare no conflicts of interest.

## Acknowledgments

We are thankful for the financial support given by the Project from Shenzhen Science and Technology Innovation Commission (JCYJ20210324105210028) and the Key Special Project for Introduced Talents Team of Southern Marine Science and Engineering Guangdong Laboratory (Guangzhou) (GML2019ZD0210 and K19313901).