Three-Dimensional Evaluation of Sand Particle Fracture Using Discrete-Element Method and Synchrotron Microcomputed Tomography Images

The fracture of sand particles plays a significant role in determining the plastic bulk volume changes of granular materials under different loading conditions. Particle breakage has been investigated using experimental, analytical, and numerical approaches at different scales ranging from a single particle to laboratory-size specimens (Altuhafi and Coop 2011; Bolton et al. 2008; Cavarretta et al. 2017; McDowell et al. 1996; Nakata et al. 1999; Zhao et al. 2015). Discrete-element modeling (DEM) has been widely used to better understand the influence of particle fracture on the behavior of granular materials. In DEM, particle fracture can be modeled using two common approaches. The first approach replaces a larger sphere with a group of smaller spheres with different diameters when a predefined fracture criterion is met (Lobo-Guerrero and Vallejo 2005; Lobo-Guerrero et al. 2006; McDowell and de Bono 2015). This approach requires selection of a breakage criterion governed by a measured or assumed characteristic particle tensile stress. The second approach simulates a sand particle as an agglomerate of small spheres bonded together using an appropriate contact model. Examples of contact models include the simple contact bond model (McDowell and Harireche 2002; Robertson and Bolton 2001), the parallel bond model (BPM) (Bobet et al. 2009; Hanley et al. 2011; Potyondy and Cundall 2004; Tomac and Gutierrez 2017), and the flat-joint contact bond model (FJM) (Potyondy 2012, 2015; Wu and Xu 2016). To better capture the actual particle morphology and microstructure of sand particles, which have been shown experimentally and computationally to play a significant role in determining the macroscopic properties of granular materials (Alshibli and Cil 2018; Andrade et al. 2012; Cho et al. 2006; Ma et al. 2017), spheres in the second approach can be replaced with bonded polyhedral blocks (Galindo-Torres et al. 2012; Gao 2013; Nicksiar and Martin 2014). The bonded block model (BBM) has been widely used in recent years to investigate fracture and fragmentation processes in rock (Mohammadnejad et al. 2018; Turichshev and Hadjigeorgiou 2017).

This paper used the BBM approach to investigate the fracture behavior of sand. Rigid blocks were used to generate zero-porosity agglomerates for a better representation of the actual morphology and microstructure of sand particles. The rigid blocks were bonded together at their contacts using the linear-parallel bond model that can transmit both moment and force (Potyondy and Cundall 2004). Each simulated particle matched the actual three-dimensional (3D) shape of sand particles acquired from high-resolution 3D synchrotron microcomputed tomography (SMT). Model parameters for the current study were calibrated based on experimental data that were acquired from unconfined compression experiments conducted on synthetic silica cubes. The model was used to compare load-displacement and fracture behavior of sand obtained from one-dimensional (1D) compression tests conducted on single sand particles and specimens composed of three sand particles.

Unconfined 1D compression experiments were conducted on three 1-mm synthetic single crystal silica cubes (MaTeck, Juelich, Germany) and two particles of ASTM 20-30 Ottawa sand that were loaded separately. ASTM 20-30 is a natural silica sand that has rounded particles and a grain size between US sieves #20 (0.841 mm) and #30 (0.595 mm). In these experiments, each sand particle or cube was placed between two loading plates and the top plate moved downward at a constant displacement rate of $0.1\mathrm{mm}/\min$ ($1.7\times {10}^{-3}{\mathrm{s}}^{-1}$ global strain rate). The SMT scans were acquired at Beamline 13D, Advanced Photon Source (APS), Argonne National Laboratory (ANL), Lemont, Illinois. Each sand particle or cube was scanned twice: before applying the load, and after fracture was observed. The acquired images had a spatial resolution of $1.98\mathrm{\mu m}/\mathrm{voxel}$. More information about the experimental setup was given by Druckrey and Alshibli (2016). Additionally, three 1D compression experiments were conducted on sand columns consisting of three particles. Sand particles were deposited into an acrylic cylindrical mold with an inner diameter of 1 mm and were compressed at a constant displacement rate of $0.2\mathrm{mm}/\min$ ($7.2\times {10}^{-4}{\mathrm{s}}^{-1}$ global strain rate). All experiments presented in this study were performed under a strain rate of less than $0.01/\mathrm{s}$, which is within the quasi-static loading condition for dry sand (Song et al. 2009). Multiple in situ SMT images were acquired at an energy level of 23 keV to produce images with a spatial resolution of $4.95\mathrm{\mu m}/\mathrm{voxel}$. Each sand column was scanned multiple times, including in the initial state before loading, during loading, and a final scan after the onset of fracture. More information about this setup was given by Cil and Alshibli (2012).

Acquired SMT scans were processed, enhanced, and segmented using Avizo 9.7 software. First, an anisotropic diffusion filter was applied to the grayscale images to reduce noise and enhance edges contrast. Then, filtered images were subtracted from unfiltered images to verify that filtering did not affect the spatial position of particles’ edges. Filtered images then were binarized using the interactive thresholding module in which user-defined values of image intensity ranges were input to binarize the images. Sand voxels were assigned a value of 1 and air voxels were assigned a value of 0. Particles then were separated using the Separate Objects module in Avizo software, and each particle was assigned a distinct numerical label. The Separate Objects module is a combination of watershed, distance transform, and numerical reconstruction algorithms that accurately remove small areas of contact between particles. More details about segmentation procedure were given by Druckrey et al. (2016). High-resolution 3D surfaces were constructed using the marching cube algorithm (Lorensen and Cline 1987) that produces a triangular approximation of the interface by computing isosurfaces from discrete data.

DEM simulations were conducted using the particle flow code PFC3D 6.0 (Itasca Consulting Group 2018), which can model polyhedral objects with triangular facets (rigid blocks). In PFC3D, the full inertia tensor of rigid blocks is used to update the rotational equations of motion. The overlap state of rigid blocks is determined using the Gilbert–Johnson–Keerthi (GJK) algorithm (Gilbert et al. 1988) which adopts the concept of the Minkowski difference to efficiently detect the overlap status of convex shapes (Itasca Consulting Group 2018). A set of rigid blocks can be constructed to fill a volume or set of volumes with a specified minimum and maximum edge length. To simulate a sand particle, rigid blocks were constructed by meshing the 3D surfaces extracted from SMT images. The fracture behavior of the sand was introduced by cementing the rigid blocks together at their contacts using the linear-parallel bond model, which can transmit both moment and force (Potyondy and Cundall 2004). This results in a crushable agglomerate of rigid blocks that closely matches the shape of an actual sand particle.

In the linear-parallel bond model, two interfaces work in parallel: a linear-elastic frictional interface, and a finite-size linear-elastic bonded interface that carries force and moment (Itasca Consulting Group 2018). Slip is introduced in the linear-elastic interface by imposing a Coulomb limit on the shear force. Because the two interfaces work in parallel, the existence of a parallel bond does not prevent slip. In PFC3D, the following parameters were used to define a parallel bond: parallel bond normal strength (${\sigma }^n$), parallel bond shear strength ($\tau$), parallel bond normal stiffness ($\overline{k^n}$), parallel bond shear stiffness ($\overline{k^s}$), contact normal stiffness ($k^n$), contact shear stiffness ($k^s$), parallel bond radius multiplier ($\lambda$), and interparticle friction coefficient ($\mu$). The bond can be pictured as a cylindrical cementing glue composed of a set of elastic springs distributed over a circular cross section on the contact plane with a beam length ($\overline{L}$) approaching zero (Fig. 1) (Potyondy and Cundall 2004). The bond breaks when tensile or shear stress exceeds the specified strength parameters. The size of the bond is controlled through $\lambda$ bywhere $\overline{R}$ = bond radius; and $R^{*}$ = equivalent radius of a circle with an area equal to the contact area. The contact area is calculated as the overlapped volume between the blocks in contact divided by the penetration depth. The penetration depth is the minimum distance that the blocks in contact must be displaced along the contact normal to no longer penetrate (Itasca Consulting Group 2018). In this study, the parallel bond effective modulus ($\overline{E}$) was used to set the value of $\overline{k^n}$ aswhere $R_1$ and $R_2$ = radii of minimal enclosing spheres of blocks in contact (Fig. 1).

Results from unconfined 1D compression tests on 1-mm synthetic silica cubes composed of a single silica crystal were used to calibrate the linear-parallel model parameters. The cube was selected due to its simple geometry which eliminates the effects particle shape and morphology at contacts. A total of 6,600 rigid blocks with a maximum edge length of 100 μm were constructed within a $(\mathrm{1,000}\mathrm{\mu m}{)}^3$ cubical volume to simulate the silica cube [Fig. 2(a)]. The system was cycled under gravity to ensure that the rigid blocks were in contact with the bottom loading plate and were permitted to reach equilibrium. Loading plates were modeled as rigid walls, and the selected loading rate was slow enough to ensure a quasi-static loading condition in which no change in force-displacement behavior was detected with further reduction of loading plate speed. The corresponding loading rate was found to be $0.035\mathrm{m}/\mathrm{s}$.

The influence of varying the model parameters on the force-displacement behavior is presented in Fig. 3. Changing the bond strength affected the peak load without changing the stiffness [Fig. 3(a)]. Alternatively, increasing the ratio of normal to shear stiffness ($\kappa$) decreased the overall stiffness and increased the peak load [Fig. 3(b)]. Changing the bond effective modulus changed the stiffness without changing the peak load [Fig. 3(c)]. Increasing the friction coefficient and radius multiplier increased the stiffness and peak load [Figs. 3(d and e)].

Because the DEM parameters have combined effects, calibration of the model required conducting multiple simulations to closely match the force-displacement behavior observed in the experiments. A total of 17 simulations were performed to identify the appropriate calibration parameters listed in Table 1. First, the stiffness of the contacts between the loading plate and the blocks ($k$) was selected by matching the stiffness measured from unconfined 1D compression of a single-block cube with the stiffness measured experimentally. Fig. 3(f) displays the variation of the model results with different values of contact stiffness. The initial ductile experimental response was not matched because it was attributed mainly to deformation of microasperities and elastic deformation at contacts (Antonyuk et al. 2005; Cavarretta et al. 2017), which are not accounted for in the rigid-blocks DEM model. Consequently, the cube used in the parametric study simulations presented previously was used to obtain the calibrated values of the remaining parameters. The parallel bond effective modulus was chosen to be 70 GPa based on the silica cube effective modulus value reported by the manufacturer (MaTeck 2019). Parallel bond normal and shear strength were assumed to be equal to allow normal and shear bond failure. A flow chart of the calibration procedure is depicted in Fig. 4. Fragmentation at the corners and edges of the cube after the peak load were observed in both the SMT images and the DEM model. Fig. 2(a) presents images obtained from SMT and DEM simulations, in which fragmentation at the corners and edges after the peak load were observed in both the SMT images and the DEM model. Fig. 2(c) shows the force versus compressive displacement curves for the experiment and the model. Table 2 presents the difference in stiffness, peak load, and displacement at peak load measured from the experiments and simulations. The peak load obtained numerically was within 2.2% of that measured experimentally. The difference in the displacement at peak load was 10.9%, whereas the difference in the average stiffness was 0.8%.

After identifying the model parameters, simulations of unconfined 1D compression tests of two ASTM 20-30 Ottawa sand particles were conducted using the same model parameters. The SMT images representing the 3D surface of the sand particle were imported into PFC3D, and rigid blocks were constructed within that surface, similar to the procedure followed when modeling the cube. Similar to the cube results, the peak load predicted by the DEM model was very close to the experimental measurement (4.4% difference) (Table 2). However, the model did not well predict the average stiffness and displacement at peak load. The stiffness obtained from the DEM model was 85% higher than the stiffness measured experimentally. The reduction in stiffness was attributed to the effects of surface roughness and asperity interaction. Otsubo et al. (2015) used bender elements to study the effects of surface roughness on the stiffness of borosilicate ballotini spheres and reported that surface roughness caused a reduction in small-strain stiffness at low confining pressure. This effect diminishes as confining pressure increases. The surface roughness effects were not captured in the current DEM model; therefore, the model failed to predict the stiffness of unconfined 1D compression tests of single sand particles.

Additionally, three 1D confined compression experiments of columns of three ASTM 20-30 Ottawa sand particles were simulated. A higher number of particles and contact points and the introduced mold confinement increased the possibilities of fractured particles. To model these experiments, each rigid block that was in contact with another block from a different sand particle was assigned a linear model. The surrounding acrylic mold was modeled as a frictionless rigid wall positioned to match the SMT images. A linear model with the same properties as the linear group in Table 1 was assigned to these contacts. The sensitivity of the model to the number of rigid blocks used to simulate particles was investigated first by generating different surfaces with different numbers of faces (triangles on the surface). A higher number of faces yields a better representation of sand surface morphology. Each surface then was imported into PFC3D and 1D compression simulations were conducted. Fig. 5 presents the variation of peak load from the DEM simulations with different rigid block sizes, and shows that the peak load tended to approach a constant value when the block edge length was less than 60 μm. Thus, the block edge length (element size) was selected to be 60 μm (8% of sand particle length) in subsequent DEM simulations.

Although the model did not predict the stiffness from unconfined 1D compression tests of single sand particles well (in terms of stiffness and displacement at peak load), the model better predicted the behavior of 1D confined compression of three sand particle columns experiments. Similarly, Cil and Alshibli (2014) reported that calibrating the linear-parallel model parameters based on unconfined 1D compression tests of a single sand particle underestimated the stress–strain behavior of laboratory-size specimens, in which sand particles were modeled using agglomerates of bonded subspheres (BSM). This is attributed to the effects of surface roughness and asperity interaction. Otsubo and O’Sullivan (2018) conducted experimental and numerical assessments of the effects of surface roughness on small-strain stiffness and reported that at low confining pressure, the overall response of soil samples is governed by asperity interaction. Moreover, they reported that the sensitivity of small-strain stiffness to surface roughness decreased with higher confining pressure. In this study, the introduced confinement in the three-sand-particle column experiments minimized the effects of surface roughness on the stiffness and overall behavior, and these were predicted well in the DEM model. This highlights the advantage of selecting the silica cube geometry in calibrating the model instead of the actual sand particle to eliminate the effects of surface roughness and improve the model prediction with higher confinement.

To demonstrate the need to capture the actual morphology of sand particles, the three-particle 1D compression experiments were modeled twice, once with agglomerates representing the actual shape of the particles and once with equivalent spherical agglomerates. An equivalent sphere is defined here as a sphere that has the same volume as the particle. SMT and DEM images along with force-displacement curves measured experimentally and based on DEM models are depicted in Fig. 6. For the first test [Fig. 6(a)], the bottom particle was the first to fracture, both experimentally and in the BBM. This observation is supported by the model proposed by Hiramatsu and Oka (1966) for the maximum tensile stress within a sphere loaded diametrically ($\sigma$)where $F$ = load at failure; and $d_f$ = distance between loading plates at failure. If we define loading distance ($d^{*}$) as the average vertical distance between the upper and lower contacts (either particle-to-particle or particle-to-plate), Eq. (3) can be used to explain why the bottom particle was the first one to fracture: it has the smallest loading distance among the three particles [Fig. 7(a)], and thus the highest maximum tensile stress. Extensive fragmentation was observed from SMT images to take place in the middle region of the bottom particle, between the upper contact area (with the middle particle) and the lower contact area (with the loading plate). A similar pattern was observed in the BBM, in which the fragments were shown in different colors. Moreover, the load versus compressive displacement behavior measured experimentally matched that obtained from real-shape agglomerate very closely [Fig. 6(a)]. The difference in peak load was 8.2% and the difference in displacement at peak load was 4.3% (Table 3). On the other hand, the DEM model using equivalent spheres exhibited a higher peak load than the experiment and very different fracture patterns. Fig. 6(a) also presents an image of a DEM simulation using BSM reported by Cil and Alshibli (2015). Fig. 6(a) shows how rigid-block agglomerates better match the shape, fracture pattern, and fragmentation of sand particles than do agglomerates composed of a large number of subspheres. Additionally, rigid blocks better match the surface contact between the agglomerates and loading plate, which was shown by Cil and Alshibli (2015) to greatly affect the load corresponding to failure.

In the second test [Fig. 6(b)], the top particle was the first particle to fracture experimentally, with multiple major cracks between the upper and lower contact areas. In the BBM, the top particle also fractured with a similar fracture pattern. The peak load in the BBM was within 1.5% difference of that in the experiment. In Fig. 7(b), the top particle has the lowest number of total contact points. The low number of contact points resulted in a higher concentration of stresses, which caused the particle to fracture even though it had the highest diametric loading distance among the three particles. Similar to Test 1, equivalent spheres had a higher peak load than the real-shape agglomerates and experimental measurements. The top particle was the first to fracture in the equivalent spherical agglomerates, but the fracture pattern in the model could not be compared with the experimental one due to the different morphology. On the other hand, BSM predicted a higher peak load and displacement, and the first particle to fracture was the middle particle, which deviates from the experimental observations and the BBM results.

For the third test, the SMT images showed that the middle particle penetrated the surrounding acrylic mold wall [Fig. 6(c)]. This fixity introduced to that particle made it susceptible to higher stresses, causing it to fracture first. This experimental flaw was not observed in the other two experiments. The fixity was introduced in the DEM models by constraining the blocks in the regions that were shown to penetrate the surrounding wall in the SMT images. The first particle to fracture was the middle one in both the experiment and the BBM, the bottom one in BSM, and the top one in the equivalent spherical BBM. Both BSM and equivalent spherical BBM underestimated the peak load and the displacement at peak load. However, the BBM predicted the peak load well, with a difference of only 4.5% from the experiment. Moreover, the displacement at peak measured from the BBM was within 16.2% of the displacement obtained from the experiment.

Breakage energy refers to the energy dissipation caused by the creation of new surfaces during fracture. Breakage energy can be calculated as the ratio of the incremental input energy to the incremental change in surface area obtained from in situ SMT images (Landis et al. 2003; Zhao et al. 2015). The change in surface area of particles and fragments was measured directly from SMT images, whereas input energy usually is assumed to be equal to the work calculated from the force-displacement curve. However, the input energy dissipates not only through breakage but also through other mechanisms such as friction and slip between particles and loading plates. DEM is a useful tool that can be used to track and monitor energy dissipation through such mechanisms, which cannot be measured experimentally (Afshar et al. 2017; Wang et al. 2012). In this study, breakage energy was calculated for the three-particle DEM models and compared with the breakage energy measured experimentally. To obtain an accurate estimation of the breakage energy, slip dissipation was subtracted from the work done by the load; the same approach was followed by Afshar et al. (2017). Energy dissipating through slippage was found to range between 25% and 35% of the total work done by the external load. This is significantly lower than the 70% reported by Afshar et al. (2017) for 1D compression of an assembly of agglomerates made of spheres. In the DEM models, the change in surface area due to fracturing was measured by tracking the broken bonds between the rigid blocks. When a bond breaks, the connected rigid blocks are separated and new surfaces are generated.

Fig. 8 shows the change of the net work of load with the cumulative change in fracture surface area for the three-particle tests measured experimentally and using DEM. In Fig. 8, the energy dissipated due to slip was obtained from the BBM and then subtracted from the total input energy for both experiments and the BBM to find the net work of load. Values of breakage energy were calculated as the slope of each curve in Fig. 8 and are reported in Table 4. Experimental breakage energy ranged from 90.3 to $605.6\mathrm{N}/\mathrm{m}$. Breakage energies measured from the BBM matched well with breakage energies measured experimentally (within 15%). Table 4 presents the total number of block contacts with another block from a different sand particle for each test. For the first and second test, the energy required to initiate fracture in the sand particles increased as the number of contact points increased. Therefore, more energy was required to fragment particles with a higher number of contact points. Hence, under the same conditions of input energy (or compressive stress), particles with a lower number of contacts are expected to fracture easier. The relatively higher value of breakage energy for the third test, however, mainly is attributed to the experimental flaw mentioned previously, in which the middle particle penetrated the surrounding acrylic mold, which was modeled as fixity in the BBM. Therefore, a major part of the input energy was dissipated through this process.

The BBM approach was used to investigate the fracture behavior of sand. Each sand particle was modeled as an agglomerate of rigid blocks bonded together at their contacts using the linear-parallel contact model. This approach provides an improvement in capturing the actual shape of sand particles by utilizing high-resolution imaging techniques such as SMT. Agglomerate micromodels were calibrated based on results from unconfined 1D compression of a single silica cube. The selection of a silica cube for calibrating the model in lieu of the actual sand particle eliminated the effects of surface roughness and improved the model prediction with a higher confinement. The results demonstrated how ignoring morphology by modeling sand particles as idealized spheres fails to simulate the fracture behavior of silica sand. Investigation of particle fracture in the three-particle sand specimens revealed that contact loading condition, number of contact points, and particle geometry interaction have significant effects on fracture pattern of particles. An increase in number of contacts increases the stress required to fracture a particle. For particles with a similar number of contact points, the first particles to fracture are those with the smallest loading distance. Energy dissipating through slippage was found to range from 25% to 35% of the total work input by the external load, and thus should be considered when estimating breakage energy. Breakage energy measured from the DEM models matched well with the breakage energy measured based on SMT images. Breakage energy results suggest that more fragmentation is expected for an assembly of particles with a smaller number of contact points under the same state of compressive stress. Finally, the BBM approach was proven to successfully predict the load-displacement behavior, fractured particles, fracture pattern, and energy released due to crack initiation by direct comparison with SMT images.

This material is partially funded by the US National Science Foundation (NSF) under Grant No. CMMI-1362510. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the NSF. The authors thank Dr. Mehmet Cil for conducting sand-column experiments, and Dr. Andrew Druckrey for conducting cube-fracture experiments. The SMT images were collected using the X-ray Operations and Research Beamline Station 13-BMD at Argonne Photon Source (APS), Argonne National Laboratory. The authors thank Dr. Mark Rivers of APS for help in performing the SMT scans. They also acknowledge the support of GeoSoilEnviroCARS (Sector 13), which is supported by the National Science Foundation, Earth Sciences (EAR-1128799), and the US Department of Energy (DOE), Geosciences (DE-FG02-94ER14466). Use of the Advanced Photon Source, an Office of Science User Facility operated for the DOE Office of Science by Argonne National Laboratory, was supported by DOE under Contract No. DE-AC02-06CH11357. The authors also thank the anonymous reviewers who contributed comments and suggestions to improve this paper.