This study utilized eight distinct agglomerate shapes as research subjects. The first one, denoted as CS, was obtained from a two-dimensional projection of actual coral sand particles. The process for acquiring these projections is illustrated in Figure 1. First, photographs of the coral sand particles were taken using a smartphone, which proved to be adequate for this purpose. It was crucial to ensure that the particles did not overlap in the images to simplify the subsequent image processing. Next, the images were transformed into two-dimensional projections of the coral sand particles, and we utilized image binarization via the ImageJ program (Schneider et al., 2012) to complete this step. Following the image binarization, additional processing steps were performed, including noise removal, hole filling, and boundary extraction. Finally, the geometry files were generated using the particle boundary coordinates, and these files can be accessed from the corresponding author upon reasonable request. The other seven were ideal geometric shapes, as illustrated in Figure 1. C1 is the circle, and C2 and C3 are the ellipses with aspect ratios of 0.8 and 0.5, respectively. S1 and S2 are square-like and square, respectively. T1 and T2 are the Reuleaux triangle and equilateral triangle, respectively. The model parameters of the seven agglomerates with ideal geometric shapes were the same as those of actual coral sand (CS), and they were employed to analyze the influence of the initial particle shape.

Furthermore, the bonded agglomerate method (Cheng et al., 2003) was employed to generate crushable agglomerates in this study. In this method, crushable agglomerates are composed of many balls bonded together by breakable bonds. The crushable agglomerates were made by generating elementary balls in geometric shapes, and hexagonal close packing was adopted to avoid overlapping the elementary balls. Considering the simulation accuracy and computational efficiency, the radius of the elementary ball was set to 0.2 mm.

Three particle size parameters (i.e., maximum Feret diameter F _max, minimum Feret diameter F _min, and equivalent diameter D _e ) and two shape parameters (i.e., aspect ratio AR and sphericity S) proposed by Altuhafi et al. (2013) were employed to describe the grain characteristics in this study. Tracking the evolution of particle size and shape is one of the most difficult tasks in actual tests. However, this can be easily achieved in DEM simulations. The acquisition methods used to determine the particle size and shape distributions before and after breakage in the simulation are as follows.

To obtain the particle size and shape distributions, one must first identify the fragments. In this study, Warshell’s algorithm (Warshall, 1962) was applied to identify the fragments, the elementary particles were deemed vertices, parallel bonds were deemed edges, and each intact agglomerate was deemed a contacted graph. If an agglomerate splits into several fragments, then each fragment is deemed a maximal connected subgraph. The specific stages of identifying fragments by Warshell’s algorithm can be found in previous studies (Laufer, 2015; Liu et al., 2017). After the fragments were identified, the boundary balls of each fragment needed to be identified. Due to the arrangement of the elementary balls as hexagonal close packs, balls with fewer than six parallel bonds are regarded as boundary balls, as shown in Figure 2. Finally, the geometrical information, including the area, perimeter, particle size and shape of each fragment, can be obtained from the coordinates of the boundary balls. Arranging the boundary balls clockwise, area A can be calculated by Eq. 1. A = 1 2 ∑ i = 1 n − 1 x i y i + 1 − x i + 1 y i + 1 2 x n y 1 − x 1 y n (1)

Herein, n is the number of boundary balls, and x and y are their coordinates. The perimeter C is simply noted as 0.02n. The maximum Feret diameter F _max is the maximum distance between boundary balls. From the direction of F _max, the plane is divided into 12 directions at 15° intervals. The maximum projection of the distance between the boundary balls in this direction is considered the diameter in this direction, of which the minimum value is the minimum Feret diameter F _min. The equivalent diameter D _e is 4 A / π . The aspect ratio AR is F min / F max . The sphericity S is π D e / C . Furthermore, the particle size and shape distributions can be determined by the area of each fragment. Figure 2 shows a comparison of the initial particle size and shape distributions obtained by two-dimensional image analysis (Zeng and Liu, 2023a; Zeng and Liu, 2023b) and simulation. The particle size and shape of the crushable agglomerates in the simulation are close to those of actual coral sand particles. In addition, the equivalent diameter distributions of the seven crushable agglomerates with ideal geometric shapes are the same as those of the CS, and their particle size and shape information can be found in Table 1. F _max50, F _min50, D _e50 , AR ₅₀ and S ₅₀ are the F _max, F _min, D _e , AR, and S, respectively, corresponding to 50% content. It should be noted that when the equivalent diameter is similar, the more regular the particle is, the smaller the maximum Feret diameter, and the larger the minimum Feret diameter.

There are two methods of specimen generation in this study, both of which ignore gravity. Method A corresponds to the laboratory tests and is shown in Figure 3A. In method A, the specimens are generated by four-layer compaction to simulate the actual layered filling. The crushable agglomerates were generated in random directions and compressed by a frictionless wall to the target position. Method B is used to generate the anisotropic specimens and is shown in Figure 3B. To increase the anisotropy of the specimens, the specimens are generated by sixteen-layer compaction in Method B. The crushable agglomerates were generated at a preset angle θ and compressed by a frictionless and inclined wall toward the target position. After half of the specimen was generated, it was removed from the mold to make space for the other half to be generated. After all the crushable agglomerates were created, the plants were returned to their original position. The initial width and height of the specimens are the same as those of the actual specimens and are 101 mm and 200 mm, respectively. To exclude specimen size effects and minimize stress nonuniformities inside a test assembly, the ratio of specimen size to maximum particle size should be greater than 5 (Jamiolkowski et al., 2004; Nie et al., 2020). In this study, the ratio is approximately 10. The number of crushable agglomerates was determined by the initial void ratio n _2D of the specimen and was 673 in this study. The total number of elementary balls was 76,636, and the number of elementary balls for a single crushable agglomerate was between 110 and 120. The initial two-dimensional porosity n _2D represents the porosity between the agglomerates, excluding the pores within the agglomerates. The results can be converted from the actual three-dimensional porosity n _3D (Hoomans et al., 1996) by Eq. 2. The initial two-dimensional porosity was 0.464 for all the specimens in this study. n 3 D = 1 − 2 π 3 1 − n 2 D 3 / 2 (2)

After the specimen was created, it was isotropically compressed by four frictionless walls to a given confining pressure, and the confining pressure was set to 100, 300, 600 and 1000 kPa, consistent with the laboratory tests. A servo algorithm was employed to ensure that this process was stable and controllable. In this stage, the parallel bond strength was set to a relatively high value of 2×10¹² N/m² to ensure that no bond breakage occurred during isotropic compression. A similar approach was adopted in previous studies (Cheng et al., 2003; Liu et al., 2017). The specimens were considered in equilibrium when the stress tolerance was less than 0.1% and the ratio of the average static unbalanced force to the average contact force was less than 1×10⁻⁴. The specimens after isotropic compression are shown in Figure 3A, B. Notably, the maximum and minimum void ratios, as determined by the method proposed by Yang et al. (2012), varied among the eight different shapes of agglomerate samples, resulting in different initial relative densities for each shape. However, after isotropic compression, their relative densities became similar, ranging from 42% to 46%, due to differences in compressibility. Figure 4 shows the agglomerate orientation distribution of the CS after isotropic compression under a 300 kPa confining pressure, with similar results for other confining pressures. For the specimen created by Method A, the number of agglomerates along the horizontal direction is slightly greater than that along the vertical direction because the direction of compression is vertical. For the specimen created by method B, the inherent anisotropy of the specimens is significant, and most of the agglomerate orientations are close to the preset angle.

After isotropic compression, the normal and shear parallel bond strengths between the elementary balls are reduced to predetermined values. To ensure that the confining pressure on the left and right sides of the specimens was nearly constant during shearing, the servo control on the left and right walls remained unchanged. Subsequently, the servo control on the bottom and top walls was removed, and the specimens were sheared by moving the bottom and top walls toward each other at a constant velocity of 0.1 m/s. In addition, the index I _uf proposed by Ng (Ng, 2006) was employed to estimate whether the specimens were in a quasistatic condition during simulation. The results show that the value of I _uf is less than 1.0%, and the simulation in this study can be considered quasistatic (Fu et al., 2019).

It is noteworthy that the rigid boundaries were adopted in this study. Previous studies have indicated that numerical specimens with flexible boundaries are more representative because the boundary effect is minimized. However, the flexible boundaries were unstable and easily distorted at high confining pressures due to the high breakability and compressibility of the coral sand. Furthermore, we found that the stress‒strain relationship obtained using rigid and flexible boundaries are close to each other at low confining pressures. Our study’s objective was also not to simulate the mechanical response of coral sand under triaxial stress states. Therefore, we believe that the rigid boundary is available.

The variability of crushable agglomerates was also investigated by single-particle compression tests. The procedure of single-particle compression simulation was as follows: first, the long axis of the geometric shape was rotated to the horizontal orientation to ensure that the loading conditions were consistent with those of the laboratory tests, and the elementary balls were generated within the geometric shape. Then, the top and bottom walls were generated as rigid compression platens. Finally, compression was performed by moving the top platen downward at a fixed velocity of 0.1 m/s until particle breakage occurred.

Figure 5A shows the survival distribution curve of single-particle crushing from simulations and laboratory tests, and the fracture strengths of these agglomerates were analyzed using Weibull’s distribution (Nakata et al., 1999) by Eq. 3. P s = ⁡ exp − σ σ 0 m (3)

Herein, P _s is the survival probability of 50 agglomerates and can be calculated by Eq. 4. P s = 1 − i N + 1 (4)

Here, i is the rank position of an agglomerate when sorted in increasing order of equivalent stress, and N is the number of agglomerates, which is 50 in this study. σ is the induced tensile stress, σ ₀ is the characteristic tensile stress at which 37% of the tested particles survive, and m is the Weibull modulus. The Weibull moduli of the simulation and laboratory test are similar. This means that the variability of crushable agglomerates is close to that of actual coral sand. Furthermore, the particle strength was mainly controlled by the normal and shear parallel bond strengths, and increased with increasing the bond strengths. However, the bond strengths could not be determined by the single-particle compression tests due to the difference between the 2D DEM simulation and 3D actual test. Due to the objective of this study being to analyze the interaction between particle shape and particle breakage, the bond strengths of this study were adjusted until the particle size and shape distributions after the simulation were closed to that after the laboratory tests, as shown in Figure 5B. In this situation, the characteristic tensile stress σ ₀ of these breakable agglomerates was 61.52 MPa, which was approximately 4 times to that of actual coral particles (15.07 MPa).

In addition, the other parameters, including the friction coefficient, normal and shear parallel bond stiffness, normal and shear stiffness of balls, were determined by the triaxial consolidated-drained tests. In previous studies, the parameters of most biaxial test models were determined by referencing classical studies (Zhou et al., 2017; Fu et al., 2019; Lü et al., 2019; Zhang et al., 2020; Bisht and Das, 2021; Wu et al., 2021). However, since this study focused on coral sand, these classical parameters were clearly inapplicable, necessitating the redetermination of model parameters. The biaxial test is difficult to perform in the laboratory; therefore, the triaxial test is commonly employed to determine the parameters of biaxial test models (Alaei and Mahboubi, 2012; Liu et al., 2017; Zhou et al., 2023), as its stress state closely resembles that of the biaxial test among all laboratory tests. In addition, our study’s objective was not to simulate the mechanical response of coral sand under triaxial stress states. Instead, we aimed to derive appropriate parameters for analyzing the interaction between particle shape and particle breakage. For this purpose, laboratory triaxial consolidated-drained tests (Zeng and Liu, 2023c) were taken as the benchmark. A comparison between the DEM simulation and the laboratory test results is shown in Figure 5C, and the model parameters of the crushable agglomerates and walls used in this study can be found in Table 2.

These findings indicate that these model parameters can be used to analyze the interaction between the particle shape and particle breakage. In addition, the simulation schemes are given in Table 3, and the subsequent analysis is based on these simulation results. Notably, the eight different shapes of agglomerates exhibited the same parameters, as presented in Table 2. The variation in the normal and shear parallel bond strengths was used to analyze the influence of particle breakage on the macro- and micromechanical behavior of coral sand.

Figures 6A, B show the deviator stress‒axial strain‒volumetric strain relationships of different particle shapes under 300 kPa confining pressure for both the unbreakable and breakable particulate assemblies. The relative density before shearing varied between 43% and 45%, and the effect of the relative density was basically negligible in this analysis. Figure 6A shows that as the aspect ratio AR and sphericity S decrease, the peak deviator stress and dilation of unbreakable agglomerates increase. This phenomenon is consistent with the findings of previous studies (Lashkari et al., 2020; Nie et al., 2020; Wu et al., 2021; Xu et al., 2021). Wu et al. (2021) and Nie et al. (2020) reported that a lower particle roundness results in a higher deviatoric stress and greater dilation during volumetric change. Nguyen et al. (2020, 2021) found that the peak drained deviatoric stress and the minimum dilatancy increased with decreasing sphericity and roundness. Nie et al. (2022) reported that with the increasing of particle overall regularity, the compressibility of the sands slightly changes, while the shear strength, stress dilatancy decrease. Figure 6B shows that the peak deviator stress and dilation of breakable agglomerates are lower than those of unbreakable agglomerates, and a similar phenomenon can be found in Figures 8A, B. This means that particle breakage leads to a reduction in strength and an increase in compressibility. In addition, Figure 6B also shows that the peak deviator stresses of breakable agglomerates with different particle shapes are similar, and the irregular agglomerates have greater compressibility than the smooth agglomerates. Furthermore, Figure 6C shows the internal friction angle φ of unbreakable agglomerates and breakable agglomerates. The internal friction angle φ was calculated by sinφ = q _f /(q _f +2σ ₃), where q _f is the peak deviator stress. φ _nb , φ _b and Δφ are the internal friction angles of unbreakable and breakable agglomerates and their differences, respectively. Both φ _nb and Δφ also decrease with increasing aspect ratio AR ₅₀ and sphericity S ₅₀ , and the influence of sphericity is greater. A similar behavior was also reported in previous studies, where the internal friction angles of unbreakable agglomerates decreased with increasing roundness (Nie et al., 2020) and aspect ratio (Liu et al., 2021). In addition, the particle shape has a small influence on φ _nb . The reason may be that irregular agglomerates are more easily broken, and particle breakage weakens the effect of particle shape on shear strength.

Figure 7A shows the deviator stress‒axial strain–volumetric strain relationships of CS with different bond strengths under a 300 kPa confining pressure. The normal and shear parallel bond strengths are the same and are denoted B _s in this study. The bond strength can reflect the fragility of the agglomerates. The higher the bond strength is, the more difficult it is for the agglomerates to break. As shown in Figure 7A, the agglomerates exhibit evident softening and dilation phenomena under high bond strength. With decreasing bond strength, the particle breakage increases, the dilation gradually disappears, and the peak deviator stress gradually decreases. However, the bond strength has a small influence on the ultimate deviator stress. In addition, Figure 7B shows an approximately linear relationship between the internal friction angle and the bond strength in semilogarithmic coordinates when the agglomerates are breakable.

In this subsection, only agglomerates CS and C3 were considered due to their smaller aspect ratio and significant inherent anisotropy. Figures 8A, B show the deviator stress‒axial strain–volumetric strain relationships of CS with different depositional angles under 300 kPa confining pressure for both the unbreakable and breakable particulate assemblies. The inherent anisotropy (depositional angle) has a significant effect on the strength and deformation of the agglomerate CS regardless of whether it is unbreakable or breakable. The peak deviator stress of both unbreakable and breakable agglomerates decreases with increasing depositional angle. The compressibility of unbreakable agglomerates is greater at low depositional angles, but that of breakable agglomerates is greater at high depositional angles. The reason is that aggregates are more easily broken at low deposition angles. Additionally, the anisotropy phenomenon of breakable CS agglomerates is less obvious than that of unbreakable CS agglomerates. A similar trend can be found for C3. Figure 8C shows the relationship between the internal friction angle and depositional angle. The internal friction angle decreases with increasing depositional angle, and this phenomenon is consistent with the findings of previous laboratory tests (Oda, 1972; Zeng and Liu, 2023b) and DEM simulation (Guo et al., 2022). The difference between the 0° and 90° internal friction angles of unbreakable agglomerates is greater than that of breakable agglomerates. This means that particle breakage weakens the effect of inherent anisotropy (depositional angle) on shear strength.

Figure 9 shows the particle size distribution curves of agglomerates CS and C1 at different axial strains, including the maximum Feret diameter F _max, minimum Feret diameter F _min and equivalent diameter D _e . Considering space restrictions, other results are not given in this paper, but their trends are briefly explained. As shown in Figure 9, the particle size distribution curves of agglomerates CS and C1 change significantly during the simulation tests. The degree of variation increases with increasing axial strain, and the variation in the irregular CS agglomerates is greater than that in the regular C1 agglomerates. Furthermore, other results show that the degree of variation increases with increasing confining pressure but decreases with increasing bond strength and depositional angle, and the variation in regular agglomerates decreases. In addition, for the irregular CS agglomerates, the variation in F _max is greater than that in F _min. However, in contrast to that of CS, the F _max of regular C1 agglomerates varied more significantly. For irregular agglomerates, particle breakage is more likely to occur in the direction of F _max. In this situation, the minimum Feret diameter F _min of the fragments is similar to that of the original agglomerate, but the maximum Feret diameter F _max decreases significantly. However, for the regular agglomerates, the maximum and minimum Feret diameters are close. When particle breakage occurs, the particle size that differs most between the fragment and the original agglomerate is the minimum Feret diameter F _min, as shown in Figure 9.

The relative breakage B _r proposed by Hardin (1985) is used to quantify the particle breakage, and it was calculated by F _max, F _min and D _e in this study. Figure 10A shows the comparison between Br _Fmin or Br _De and Br _Fmax for eight different agglomerates, and the results for different confining pressures, axial strains, bond strengths and depositional angles were considered. Notably, there is an approximately linear relationship between the relative breakage B _r calculated by F _max, F _min, and D _e, regardless of the agglomerate shape, and the slope is related to the aspect ratio AR ₅₀ and sphericity S ₅₀ . In this study, the slopes k _mm and k _em of Br _Fmin -Br _Fmax and Br _De -Br _Fmax can be estimated by Eq. 5 and Eq. 6. k m m = 0.475 + 1.117 A R 50 − 0.258 S 50 (5) k e m = 0.485 + 0.330 A R 50 + 0.448 S 50 (6)

Herein, AR ₅₀ and S ₅₀ are the initial AR and S corresponding to a 50% content, respectively.

Figure 10B shows the relationship between Br _bn and Br _Fmax . Br _bn is the relative breakage calculated from the number of broken parallel bonds by Eq. 7. B r b n = N b r o k e n N a l l (7)

Herein, N _broken is the number of broken parallel bonds, and N _all is the number of parallel bonds without particle breakage. Br _bn and Br _Fmax are approximately equal for the CS agglomerates. In addition, Figure 10C indicates that Br _De can be estimated by Br _Fmin and Br _Fmax , and their relationship is Eq. 8. B r D e = 1 − 1 − B r F m a x 1 − B r F m i n (8)

For the simulation tests, the input energy per unit volume E from the start of the test (SOT) to the end of the test (EOT) can be expressed by Eq. 9. E = ∑ S O T E O T σ 1 − σ 3 d ε 1 + σ 3 d ε v (9) where σ ₁ and σ ₃ are the major and minor effective principal stresses, respectively, and dε ₁ and dε _v are the axial strain and volumetric strain increments, respectively. The isotropic compression process was also considered in the analyses, and a confining pressure of 10 kPa was used as the initial condition. Although the agglomerates were unbreakable during isotropic compression, the breakage caused by the reduction in bond strength before shearing can be regarded as particle breakage during isotropic compression. As shown in Figures 11A, B unique relationship can be established between the relative breakage Br _De and the input energy E, regardless of the axial strain and confining pressure, but this relationship is significantly influenced by the agglomerate shape and depositional angle. Under the same input energy E, the relative breakage Br _De increases with decreasing depositional angle, and the agglomerates with small aspect ratios AR and small sphericity S are more easily broken.

Figure 12A shows the particle shape distribution curves of agglomerates CS and C1 at different axial strains, including the aspect ratio AR and sphericity S. Considering space restrictions, other results are not given in this paper, and the trends are similar. As shown in Figure 12A, the particle shape distribution curves of agglomerates CS and C1 also change significantly during the simulation tests. However, the trends of agglomerates CS and C1 are obviously different. For the irregular CS agglomerates, with increasing axial strain, the aspect ratio AR and sphericity S increase, and the particles become more regular. However, for the regular C1 agglomerate, the particles become more irregular with increasing axial strain.

The shape factors Sr _AR and Sr _S were employed to further analyze the evolution of the particle shape. Particles with aspect ratios less than 0.2 or sphericities less than 0.5 are basically nonexistent. Therefore, the shape factor Sr _AR was defined as the ratio of the area between the shape distribution curve and the vertical line with an AR of 0.2 to the area between the vertical lines with an AR of 0.2 and 1. The shape factor Sr _S was defined as the ratio of the area between the shape distribution curve and the vertical line with an S of 0.5 to the area between the vertical lines with an S of 0.5 and 1. A larger shape factor represents a more regular particle shape. Figures 12B,C show the evolution of Sr _AR and Sr _S during simulation, respectively. The evolution of the particle shape was controlled by the particle breakage, regardless of the axial strain, confining pressure, bond strength and depositional angle, and the trends were determined by the initial shape factors Sr _AR and Sr _S . The shape factor Sr _AR gradually tends to 0.7 with increasing particle breakage. The shape factor Sr _AR increases and decreases when the initial shape factor Sr _AR is less than or greater than 0.7, respectively. When the initial shape factor Sr _AR is close to 0.7, the shape factor Sr _AR changes slightly. The shape factor Sr _S decreases first and then increases when the initial shape factor Sr _S is greater than 0.75, and it has a significant increasing tendency when the initial shape factor Sr _S is less than 0.75.

Based on the analysis above, the interaction between the particle shape and particle breakage is highly significant. The particle shape affects the difficulty of particle breakage, and particle breakage controls the evolution of the particle shape. In addition, they have combined effects on the mechanical behavior of coral sand.