Effect of particle morphology and contacts on particle breakage in a granular assembly studied using X‑ray tomography

The macroscopic response of a geomaterial is entirely determined by changes at the particle scale. It has been established that particle crushing is affected by particle size, shape and mineral composition and initial density; and the initiation of breakage has often been related to the onset of yielding. Encouraged by the success of X-ray tomography in revealing particle-scale mechanisms of deformation, we present our findings regarding the onset of particle breakage, deriving from our study of 3D images of a dry granular assembly undergoing crushing. We propose two bespoke image analysis algorithms, which allow us to track breakage and identify contacts prior to breakage. The combination of the two algorithms, along with the high resolution of the 3D images enables us for the first time to track breakage of individual particles, identify different breakage modes for each particle and simultaneously study the effect of particle morphology and coordination number on breakage. Three different breakage types are identified: chipping, splitting and fragmentation. We have found that particle heterogeneity and sphericity mainly contribute to fragmentation, whereas the coordination number also affects chipping. The confining stress state within the particles with high coordination number made them more resistive to fragmentation, whereas particles with low coordination number mainly undergo fragmentation. The shearing of the particles at their contact points, leads to local stress concentrations resulting into surface chipping. Finally, we discuss the relation between the initiation of breakage and yielding, showing that some breakage occurs before the point where yielding is traditionally defined.


Introduction
In the study of grain breakage phenomena in particulate matter, standard quantitative measurements are usually limited to an initial and a final sieving analysis, and to the externally-measured response of a specimen to applied loading. This very limited information is a hindrance to researchers aiming to inject accurate micro-observations within a micromacro modelling of phenomena such as breakage, comminution or volumetric collapse.
Particle crushing is affected by particle size, initial density, particle shape and mineral composition (e.g., [1][2][3][4]). In a number of studies (some presented in the following), density, specific volume and/or void ratio are used to study the evolution of breakage, often implying the effect that fabric, and specifically coordination number, has on breakage.
In a granular assembly, particles break with increasing applied macroscopic stress, and particles with lower coordination number and larger particle size are found to The work presented in this manuscript is a continuation of the doctoral work of Dr. Karatza, part of which has been presented in scientific conferences (references have been made throughout the manuscript). be more prone to breakage [5]. On the contrary, numerical analyses [5] and experimental findings supported by DEM evidence using crushable agglomerates [6] show that coarse particles survive in a granular material undergoing crushing [4]. This is because in a polydisperse assembly a decrease in coordination number and an increase in particle size are two opposing conditions. McDowel et al. [5] explained that crushing tensile strength is expected to fall with increasing particle size due to the increased likelihood of internal flaws existing in the larger particles, yet smaller particles tend to break as increased coordination number outweighs the effect of particle strength solely due to particle size (also reported by [7,8]).
Nakata et al. [1] showed that the average characteristic tensile stress sp in an assembly ( sp = [(1 + e) ∕6] 2∕3 , where is the applied stress and e the void ratio), at which a particle would fail, is independent of its size, but is a function of the void ratio of the assembly. Larger particles in an assembly, having a higher coordination number [4], get cushioned by surrounding smaller particles, making them more resistive to crushing. Therefore, finer particles are more likely to be crushed. This cushioning effect was observed numerically and experimentally by [9] and also described by [10] and [2]. Finer particles are subjected to loading similar to that of a Brazilian test performed on concrete, due to the small coordination number, resulting in tensile failure by splitting. Coarser particles have a reduced probability of such failure, as the increasing coordination number results in a confined stress state within the particles [5]. Contrarily, a recent experimental study by [11], showed that for particles with a lower coordination number (specifically 4) failure occurred by splitting, fragmenting or chipping, but for those with a higher coordination number (specifically 6) only splitting occurred.
Particle composition also affects the fracture patterns in a grain. Irregular morphology influences the stress distribution inside a particle; local stress concentrations can occur at contact points with small radii of curvature, leading to extensive fragmentation [12].
The term yielding is used to describe the post-elastic behaviour of granular media, and in oedometric compression tests it is often related to successive particle breakage [1,5].
McDowell [3] explains that yielding is a gradual process, and it is the yield stress (i.e., onset of yielding) that studies refer to as yielding, beyond which an almost linear normal compression line emerges. In oedometric compression tests, yielding is usually identified on the e − log plot as either the point of maximum curvature [3], the point where the relation of the deformation to the stress increment changes rapidly [1] (i.e., turning region with stiffness reduction) or the point defined by Casagrande [13], all of which can be determined analytically or graphically.
Coop and Lee [14] performed oedometric compression tests on loose specimens of three different sands and showed that the stiffness reduced due to particle crushing and attributed yield to the onset of crushing. They found that the yield stress depends on the strength of the soil particles and the relative density of the specimen. Others have also attributed yielding to the onset of particle breakage, defined as soon as the stiffness starts to change [1,4,5,[15][16][17][18][19][20][21]. However, McDowell [3] states that "at yield stress the particle size distribution will already be different from the initial distribution...", implying that significant breakage will occur before yielding, in order to have evident changes in the particle size distribution. Finally, it is found that the yield stress will increase with decreasing particle size [1,22] and increasing initial void ratio [23] and will be approximately proportional to the tensile strength of the constituent grains [7,22,24,25].
This paper aims to make several contributions. Encouraged by the success of X-ray tomography in revealing particle rearrangements during shearing tests (e.g., [26][27][28][29]), as well as revealing some parts of breakage processes at higher mean stresses (e.g., [30][31][32][33][34][35][36]), we present our findings regarding the onset of particle breakage, deriving from our study of 3D images of a dry granular assembly undergoing crushing. Specifically, strain controlled oedometric compression tests were carried out on an assembly of zeolite granules, while performing X-ray tomography, focusing on the initial stages of particle breakage. We propose two bespoke image analysis algorithms to track breakage, measure a number of particle properties (size, shape, composition), as well as the number of contacts per particle prior to breakage. This work aims to provide a deeper understanding on the effect that coordination number and particle morphology have on different breakage modes, and on their contribution to the onset of breakage. Finally, we discuss the relation between the initiation of breakage and yielding, as it has been observed in this particular set of experiments.

Experimental campaign and macroscopic results
Industrially manufactured zeolite granules are chosen for this study due to their high roundness and sphericity. The particles have a density of 2.18 g cm −3 and an estimated crushing strength of 15 N (based on the supplier's specifications). The mean diameter ( D 50 ) of the tested sample is 1.36 mm , the minimum diameter ( D min ) is 1.09 mm and the maximum diameter ( D max ) is 1.50 mm . The grading is very uniform (coefficient of uniformity C u = D 60 ∕D 10 = 1.07 ) allowing all intact particles to be imaged with the same detail, when X-ray computed tomography (XCT) is performed. The particle assembly has an average sphericity of 0.97, which is measured directly from the 3D images, as a function of particle volume V and surface area A, as follows: The sphericity shows how well a particle can fit into a sphere, which means that spheres have equal to one and elongated ellipsoids have closer to zero, despite both being highly rounded. The zeolite granules are tested under strain-controlled oedometric compression. In this paper results from one test are presented (test OCS1) in which XCT is performed; however, a number of repeat tests were run (see [37]) to ensure the repeatability and reproducibility of the experimental procedure, of which two had XCT performed (OCB4 and OCS2). In the naming of the different tests, OC refers to oedometric compression and B or S to the size of the oedometric cell (B referring to a larger cell, with 45 mm diameter and 15 mm height). For test OCS1, the height of the specimens is 15 mm , chosen to satisfy [38] standards for one-dimensional consolidation test, which relates the size of particles to the height of specimen, and the slenderness ratio equal to 1. The specimens were created to have the same dense initial configuration (porosity, n ≈ 40% ). The oedometric cell was specially designed for the needs of this study. It is made of PEEK (PolyEther Ether Ketone), chosen for its low X-ray absorption and low friction. Figure 1a shows a schematic of the cell and of the filling and loading. The loading direction is upward and the quasi-static loading is performed with an axial loading rate of 50 μm min −1 . Figure 1b shows the response of the tests in which XCT was performed.
The macro-mechanical behaviour of the tests, as presented on an e − log plot, follows the traditional threepart response. Initially we observe a linear stiff response, followed by a turning region after which we notice an almost linear normal compression line. At higher stresses (beyond 5 MPa), due to the excessive breakage, the void ratio has (1) = 1 3 (6V) 2 3 A become very low and so the curvature changes again. This change of curvature is explained by [3] as a change in the material making it behave more like a rock (increase in stiffness). The remaining large particles (intact) are well protected by the smaller neighbours (fragments), so the tensile stresses induced within the large particles are very low, while the smallest particles have reached the comminution limit and so additional fracture is no longer possible, increasing the stiffness of the specimen.
The tests were performed in Laboratoire 3SR (Grenoble, France) using an X-ray scanner with 100 kV of acceleration (described in [39]). The voxel size is 24.9 μm, 12.3 μm and 10.0 μm for tests OCB4, OCS1 and OCS2 respectively, providing images of high resolution (more than 100 pixels per intact particle diameter for test OCS1, results of which are presented here). Loading is performed in-situ in the scanner under displacement control, setting the strain rate to zero during scanning. After inspecting all three tests, the results regarding breakage and the macroscopic response are consistent, which is why only results for OCS1 are presented here. Figure 2 shows vertical slices and details of each load step (LS#) of test OCS1. Initially breakage is observed close to the moving boundary and it is related to the friction between the loading platen and the particles in contact with it, as well as the filling process (see Fig. 1a).
Visual inspection of Fig. 2 shows that grain breakage produces fragments that are either big enough to be clearly identified as new particles, or too small with respect to the resolution of the images and are just taken as "fines" (< 0.012 mm in diameter). In the first part of the response, before the turning region (LS1 and LS2), there is no breakage. When particles start to break (LS3, just before the turning region) all of the breakage is concentrated within two to three intact particle diameters of the bottom boundary. Right after (LS4) the number of particles breaking close to the moving boundary increases and some breakage occurs in the rest of the specimen. Until LS4, breakage occurs only to intact particles (primary breakage), whereas right after LS4 fragments also start to break (secondary breakage). Until LS7 (second change of curvature), there are hardly any "fines". After LS7, there is an increase in the amount of breakage of both intact particles and their fragments. In LS9, mainly fragments continue to break, yet after LS10 there is no indication of secondary breakage, just reduction of the void space and most of the intact particles identified in LS10 are also identified as intact in LS11, indicating the material has reached a comminution limit consistent with the findings of [3]. From the inspection of the images it can be seen that breakage occurs at the start of the turning region, instead of the point of maximum curvature on the e − log plot, as suggested in the literature. Quantitative results and further discussion are presented next. The main focus of this paper is the initiation of breakage, therefore in the following, results are presented for the first three loading stages of OCS1.

Image processing: tools and algorithms
The use of XCT, has improved the study of micro-mechanics, allowing 3D visualisation of processes and 3D full-field measurements [40]. In this work, two algorithms have been developed to track and quantify breakage, and relate it to a number of particle properties, measured directly from the 3D images. Initially, all the intact particles are identified in the initial 3D image and a number of shape and size properties are measured. An algorithm was developed to track the intact particles in the subsequent images and therefore to quantify breakage. To relate the coordination number to breakage, a second algorithm was developed to identify particle contacts. Both algorithms are described in the following.

Discrete digital image correlation
Digital Image Correlation (DIC) is an image analysis method, applied to non-contact full-field kinematics measurements of planar or non-planar surfaces (2D DIC) or volumes (3D DIC) undergoing deformation [41][42][43]. In 3D DIC, a pair of images is compared to each other, where the initial image is separated into cubic sections (correlation windows) and each of them is correlated with cubic sections of the same size within a predefined region (search window) in the second (deformed) image. In the deformed image, the correlation window with the highest normalised cross correlation of the mean greyscale value of all correlation windows within the search window, represents the new state of the cube under investigation. Based on the spatial coordinates of the undeformed and deformed correlation windows, displacement vectors are defined [44].
In this work, the "discrete" mode of TomoWarp2 [45] (dDIC) is used, where correlation windows are individually adapted to each particle and only the voxels corresponding to each particle are used in the correlation. This method is essentially the same as the one presented in [27]. By taking a total approach from the beginning of the test (initial image) tracking intact particles into the following images, broken particles are defined as particles where the correlation fails (i.e., the highest correlation coefficient is below a certain threshold). Briefly the algorithm is described in the following:

Initial grain identification
The solid phase of the initial image is separated into individual particles using a standard watershed [46]; despite the documented errors with such an approach (e.g., [47]), the highly spherical particles are not a challenge to segment and label accurately. The segmentation is performed with an accuracy of 99.8% (validated against the mass of the specimen and the D 50 of the particles, both measured independently prior to testing). A result of the segmentation is that in the initial image each particle has an individual label. The advantage of the dDIC is that the tracking of the intact particles does not require the segmentation of the deformed images, where the segmentation can be poor due to the low void ratio and small size of fragments.
Tracking particles After defining the correlation window based on each particle's shape, the centroid and the mean and SD of the greyscale value of each particle, are also determined. The correlation is performed using a pixel search followed by sub-pixel refinement by image interpolation (as described in [45]). The shape and centroid of each correlation window are used to track each particle's movement (displacement and rotation).
Identifying particles that break After tracking each particle in the next image, any grain whose sub-pixel normalised correlation coefficient is lower than a given threshold (in this case 96%) is taken as being broken (see also [36]). This threshold was carefully selected after studying the distribution of the correlation coefficients at the end of each dDIC.

Re-labelling grains in the deformed images
An algorithm is currently being developed to precisely relate fragments to intact labelled particles. Due to the small number of broken particles in LS2 and LS3, however, simply identifying fragments in the space that would be occupied by an intact particle, was enough to relate those fragments to an intact labelled particle. This required the deformed images to be segmented, so that each labelled particle (intact or fragment) in the deformed image can be re-labelled, using the label of the intact particle, which it relates to. This step is not currently a feature of the dDIC. Figure 3 summarises the key steps of the dDIC.

Contact detection
Commonly, in image processing for granular media, the image is binarised (i.e., separated into two continuous phases) and segmented (i.e., binarised grains separated from each other) and the contact area is defined as the resulting voxels from the watershed plane when the two images (binary and segmented) are subtracted from each other (e.g., [48,49]). Limitations can arise depending on the resolution of the images, the noise, the partial volume effect and the choice of the thresholding and segmentation algorithms. Due to these limitations, problems in the detection of contacts and the determination of their orientation may arise by overdetecting contacts or having an unrealistic orientation (see [50,51]).
Here a different approach is proposed, using the gradient of the greyscale images to refine the contact detection described previously. Initially, an estimate of the contact area is calculated from the subtraction of the binarised and segmented images. The greyscale image is then converted into a continuous function (f(x, y, z)) and the magnitude of the gradient of the continuous image ( (x, y, z) ) is computed, as follows: where x, y and z are the coordinates of each voxel. A threshold is defined for the greyscale image of the gradient, to separate the areas with a high gradient from the rest of the image. Voxels with a high gradient indicate a transition from a void (low greyscale value) to a particle (high greyscale value). Voxels of contacting particles have a slight (if any) change in greyscale value and so locally, at the contact area, the gradient is close to zero. This is then used as a criterion to verify and/or alter the contact areas that have been defined from the subtraction of the binary and segmented images. The thresholded image of the gradient will either show a continuous line around each particle (i.e., high gradient = no contact) or it will be interrupted at the area of an existing contact (i.e., low gradient). If therefore we can virtually move from one particle to the other, without crossing a particle boundary, then we are moving through an actual contact. The existence of particle boundary voxels in the image of the thresholded gradient, at the location of an estimated contact (binary-segmented), results in this contact being removed or reduced in size.
The proposed contact detection algorithm will be referred to as the gradient algorithm, whereas the algorithm detecting contacts by subtracting the segmented image from the binary will be referred to as the bin-segm algorithm hereinafter. Figure 4 summarises the key steps of the contact detection algorithm.

Fig. 3 Schematic of discrete Digital Image Correlation algorithm.
[Note: in the Correlation Results, the brightness and contrast have been adjusted, so that all particles with a CC value higher than 96% appear with the same dark colour and all those with a CC value lower than 80% also appear with the same light colour.]

Fig. 4 Schematic of contact detection algorithm
Validation schemes A DEM model of test OCS1 was created, using material properties calibrated based on zeolite's response to single particle compression, and contact and physical parameters to match the initial configuration of OCS1 in terms of packing, grading and number of particles (described in [31]). The distribution of the number of contacts per particle in the DEM model matches the one being measured from the images, using the gradient algorithm. Instead, the mean and the SD were overestimated, when identifying the contact areas with the bin-segm algorithm. The range of coordination numbers found in this work using the gradient algorithm is also consistent with ranges reported in the literature for dense configurations of polydisperse assemblies of spherical particles (e.g., [52]).
Additionally, fabricated images of two spheres (contacting or not) were created using the Kalisphera tool [53] to investigate the performance of both contact detection algorithms. For OCS1 the image pixel size is 12.2486 μm (res) and the smallest intact particle size is 1.09 mm ( D min ), which were respectively the resolution and sphere diameter for the fabricated spheres. A set of ten images were created, where the spheres were overlapping (particles' centroid distance equal to D min − 0.4res and D min − 0.2res ), in point contact (particles' centroid distance equal to D min ) or not in contact. In the case where the particles are not in contact, the centroids' distance was set to be one particle diameter plus 0.2 to 1.2 of the image resolution, in steps of 0.2 of the resolution, to investigate the effect of the partial volume effect (see [54,55]) on the contact detection algorithms. In addition, noise was introduced to the images, to investigate the sensitivity of the algorithm to noise, and blur was also added to measure the resolution of the algorithms (not to be confused with the resolution of the images, which refers to the pixel size). For the binarisation we used the method proposed by [56], as when it was compared to other thresholding algorithms it returned the most accurate results in terms of the volume of the spheres. It should be mentioned that the gradient algorithm however, is not affected by the thresholding method. Figure 5 shows a comparison of the two algorithms with respect to blur (left) and noise (right). In the legend the first two letters indicate the algorithm (BS for bin-segm and Gr for gradient), the following letter refers to the addition of blur (B) or noise (N) and the number shows the amount of blur or noise that has been added to the images. Both figures have the same two dashed vertical lines, the left one indicating the position of the point contact where the contact area should be 0 mm 2 and the right one showing the position of the particles when they are one pixel apart.
As it can be seen in Fig. 5 both algorithms correctly identify the non-existence of a contact when the distance between the particles is equal to or larger than one pixel. This is due to the high resolution of the fabricated images, in terms of sphere and pixel sizes. For any algorithm to detect a point contact, at least one pixel has to belong to the contact. This means that the minimum contact area for a point contact is equal to the square of the resolution ( 150 μm 2 ). While both algorithms overestimate the contact area, the measured area using the bin-segm is almost twice the one measured using the gradient algorithm. Results from the gradient algorithm show that it is more sensitive to noise but it is affected less by the blur in comparison with the bin-segm algorithm. Most importantly, the gradient algorithm successfully identified the existence of a contact for reasonable values of the blur and regardless of the noise, validating our choice of using the gradient to develop this new algorithm. Using the bin-segm algorithm for this study would significantly overestimate the coordination number in OCS1, as it is shown in the results presented in Fig. 5. It should be mentioned that errors could be introduced during the thresholding and segmentation of the fabricated images, which are not evaluated at this stage.

Particle-scale results
The combination of the dDIC and the contact detection and the high resolution of the 3D images allows us for the first time to simultaneously track breakage of individual particles  (temporally and spatially), identify different breakage modes for each particle and study the effect of particle morphology and coordination number on breakage. Figure 6 shows the evolution of broken particles that are identified in OCS1 (from LS2 to LS5) using the dDIC algorithm. Since the main focus of this paper is the initiation of breakage, and a large number of broken particles are identified after LS3, only results for LS2 and LS3 are analysed in depth hereafter. Table 1 presents a summary of measurements with respect to breakage for test OCS1. In this test 1058 intact particles are identified at the beginning of loading, 0.2% of which break in LS2 and 8% in LS3 (just at the beginning of the turning region). This clearly indicates that at the maximum point of curvature on the e − log plot (where yielding is typically identified) a significant amount of particles will have already broken. A similar response to breakage was measured in the other tests, where XCT was performed (OCS1 and OCB4). Particles are classified, manually, into three breakage modes: chipping, splitting and fragmentation. Figure 7 shows the distribution within the sample of the different breakage modes. Only four particles (5% of the broken particles) split axially, therefore the results presented for this category might be statistically insignificant. The majority of broken particles (54%) broke by fragmentation, particularly those near the bottom cell boundary (a similar boundary effect is reported by [57]).
To explore the effect of particle morphology on the breakage pattern, Fig. 8 shows histograms summarising the effect of particle heterogeneity, size and shape including all 1058 intact particles identified in LS1. The choice of graphical representation of the results was investigated extensively, as it was found to significantly affect our ability to draw general conclusions. A typical histogram, with constant bin width, would lump most measurements in very few bins. This would leave the remaining bins with few occurrences, and therefore highly susceptible to statistical error. For this reason, a histogram with equal number of occurrences (i.e., height) per bin was chosen. Figure 8 presents results in five bins (i.e., quintiles) of ∼ 212 particles each, resulting in nonconstant bin widths. The overall height of each bin would be 212, however we have removed the results of the particles that remain intact in LS3, to better visualise the effect of the investigated parameters on breakage.
Shape is studied in terms of particle sphericity (see Eq. 1), regardless of whether the particles are angular or rounded. Figure 8a shows a clear effect of the particle shape on breakage, more so than the other parameters. The appropriateness of our choice of visualisation is clearly demonstrated in this plot. It is immediately evident that the lower quintile of the particles, in terms of sphericity, exhibits the most breakage (especially fragmentation), while the upper quintile exhibits almost no breakage at all, contrary to the   Fig. 7 3D image of intact particles in LS1 that will break in LS3, categorising them based on the breakage mode findings of [11]. It should be mentioned that the upper limit of the particles' sphericity is restricted by the fact that the particles are pixelated and even a perfect sphere will not have a sphericity of 1, as expected. No correlation was found between particles of low sphericity and high heterogeneity or low density, which could justify why less spherical particles tend to break more. Therefore it is suggested that particles with lower sphericity have a smaller radius of curvature at their contact points with neighbouring particles, leading to local stress concentrations at the contact points causing fragmentation (also observed by [12]).
To study the effect of the size of the particles, the particle volume is used, shown in Fig. 8b. Note that the same study was performed using the particles' D 50 and the results were consistent, largely because of the high sphericity of the intact particles and the narrow grading of the intact assembly. Many studies have reported that the characteristic tensile stress induced within a diametrically loaded grain is inversely proportional to the particle's diameter and its flaws (e.g., [1,5,58]), leading to the conclusion that larger particles will fail more easily when compared to smaller particles of the same mineral composition and number of flaws. It is also suggested in the literature that smaller particles will have fewer flaws than larger particles; no such conclusion can be drawn from this study, however there is only a small difference in size between the smallest and largest particles in the specimen. Figure 8b suggests that particle size does not have a strong influence on breakage. There is indeed a larger number of breakage for the larger particles, but equally so for the particles of lesser volume (especially in fragmentation mode).
Particle heterogeneity is measured in terms of the greyscale value of the voxels belonging to each labelled particle. Low average particle greyscale values (i.e., low particle density) suggest the presence of high internal porosity (flaws in the form of pre-existing hairline cracks), whereas high values suggest dense inclusions (see white spots in particles in Fig. 2). The greyscale values can vary between 0 and 65535 in a 16-bit unsigned integer format image. Here these values are normalised between 0 and 1 to allow an easier interpretation of the plots showing results for the mean and SD of the greyscale values. The combination of the average and the scatter of greyscale value of broken particles (Fig. 8c, d) can provide some insight on whether and if so how particle density affects breakage. As it can be seen in Fig. 8c and d, the heterogeneity of the particles mainly affects the fragmented grains. Specifically, particles with the highest heterogeneity (i.e, high SD- Fig. 8c) are those that tend to break more, particularly those with a lower average greyscale value (Fig. 8d) indicating a lower overall particle density. The combination of the two suggests that particles with an overall lower particle density and a larger variation of internal composition arising probably from internal flaws are the ones most likely to break.
To explore the effect of the coordination number on the breakage mode, the average coordination number per Fig. 8 Effect on breakage of a particle shape; b particle size; c particle heterogeneity (SD of greyscale value ranging from 0 to 1) and d particle density (average greyscale value ranging from 0 to 1) in LS3. Note that each figure has ∼ 212 particles in each range, including the intact particles 44 Page 10 of 13 category is measured and presented in Table 2. The particles with more contacts undergo surface chipping and the ones with the fewest exhibit a breakage behaviour resembling that of a Brazilian test (i.e., axial splitting). The 3D images shown in Table 2 also contribute to the understanding of the spatial arrangement of the contacts, relating them to the various breakage planes, however, an analysis as such is not presented in this work. To better understand the effect of the coordination number on the breakage mode, results are shown in Fig. 9. Figure 9a, shows a typical histogram, including results of both intact and broken particles for each breakage mode. The small number of broken particles makes it difficult to investigate the effect of the coordination number on breakage. Therefore, the histogram in Fig. 9b was created where the data in each bin is normalised with the total number of particles each category so that the height of each bin reaches 100%, if results for the intact particles are added to each range. The results show that particles with lower coordination numbers exhibit mainly splitting and fragmentation, whereas chipping occurs to particles with higher coordination numbers. Our findings are contrary to those of [11], however, complementing them as they have only investigated particles with two, four and six contacts. Overall, there are two competing general failure modes: surface failure (chipping) or internal failure (fragmentation). In cases where particles have a high number of contacts, the internal tensile cracks initiation can be curtailed due to confinement from the multiple contacts (provided that the contacts carry strong forces; see definition of [59]), resulting in surface failure being dominant.

Conclusions
In this paper, the evolution of particle breakage is investigated during strain-controlled oedometric compression tests, performed on well-rounded, highly spherical and uniform zeolite particle assemblies. X-ray computed tomography is employed to acquire high-resolution 3D images of a number of specimens at different stages of loading. The experiments  are found to be highly repeatable and the results from one test has been analysed in depth and presented. Two bespoke algorithms are presented, which allows us for the first time to simultaneously track particles and their fragments in the initial stages of breakage and investigate how particle morphology and coordination number affect the different types of breakage. The main findings are: • Visual inspection of the images shows that while in the initial stages of breakage only intact particles break, once loading progresses both intact particles and fragments continuously break. At the end of the test with a vertical stress of 25 MPa, a number of zeolite particles still remain intact. • A significant amount of breakage is detected well before the typically defined yielding point, that is often associated with the onset of particle breakage. • Three different breakage types are identified; chipping, splitting and fragmentation. While splitting occurred to very few particles in the test, there was an almost equal contribution from the remaining breakage modes. • There is a clear boundary effect, causing the majority of breakage to initiate closer to the loading platen and gradually propagate upwards to the rest of the specimen. Interestingly, all of the particles in contact with either of the top and bottom boundaries exhibited fragmentation. • Particle heterogeneity (internal density scatter) and sphericity mainly contributed to fragmentation, whereas the coordination number also affected chipping. • The confined stress state within the particles with high coordination number made them more resistant to fragmentation. When a particle has a high coordination number, the internal tensile cracks initiation is being curtailed due to confinement from the multiple contacts, resulting in surface chipping being dominant. • Particles with low coordination number undergo fragmentation and splitting, as a result of the compressive forces imposed at their contacts with their neighbours.
Experimental quantifiable evidence is presented in this work, regarding particle-scale parameters that lead to particle breakage. Such results could lead to improved continuum models. A limitation of the current work is presented in the fact that fabric is not fully investigated, even though it is known that the magnitude and orientation of the contact forces greatly affect the way particles break. Contact forces have not yet been confidently measured directly from XCT data (see [29,[60][61][62]), which is why they are not presented in this paper. Our experimental work will be coupled with a DEM analysis of test OCS1. The DEM parameters are calibrated based on the material response and an initial configuration is created that matches the experiment in terms of packing, grading and number of particles. This coupling will allow a further investigation on how fabric and particle morphology affect particle breakage at its early stages.