DEM studies on the effect of particle breakage on the critical state behaviours of granular soils under undrained shear conditions

A number of discrete element analyses of undrained triaxial shear tests on crushable assemblies are performed using 3-dimensional particle flow code (PFC3D). The undrained shear of the samples is simulated by assigning different velocities at the boundaries to control the constant volume of the samples during the shear. Particle breakage is studied using the octahedral shear stress breakage criterion and the fragment spawning mode satisfying an Apollonian sphere packing. The microparameters of the soils are obtained by simulating the test results reported on a Nansha calcareous sand. The effects of particle breakage and drainage condition on the critical state behaviours of the sand are intensively examined. It is found that particle breakage and drainage condition do not affect the slope of the critical state p′−q line. In the ecs−(p′cs/pa)0.7 space, the critical state line translates downward and rotates clockwise when particle breakage is considered, whilst drainage condition does not affect the critical state line. The intrinsic macro–micro correlations underlying the critical state behaviours are explored by analysing the evolution of representative micromechanical indicators such as coordination number, redundancy ratio, and fabric and force anisotropies.


Introduction
The pioneering work by Roscoe et al. [65] on the critical state soil mechanics (CSSM) framework lays a footstone for the development of modern soil mechanics [3,5,66,82]. Up to date, many CSSM-based constitutive models have been successfully used in analysing soil behaviour [27,31,41,43,95]. The key concept of CSSM is that when the soil reaches the critical state, the soil deforms at a constant volume under constant resistance during shear. Many researchers have discussed the uniqueness of the critical state line (CSL) in the e-p'-q space (where e is the void ratio, p 0 is the mean effective principal stress and q is the generalized shear stress). In the earlier research, many factors affecting the uniqueness of CSL have been examined, such as soil grading [29,32,89,102], particle shape and surface asperity [10,35,60,72,90], and stress path [57,87,100]. Recently, the issues of particle breakage have attracted the attention of researchers [15,17,52,68,73,78,99], due to the frequent encounter of crushable materials in geotechnical engineering such as calcareous sand, which can crush at very low stress level, e.g. B 1 MPa [14]. Furthermore, the effect of particle breakage on the critical state behaviours of granular soils has also been investigated by many researchers using different methods [4,20,30,74].
Experimental investigation is an important approach to studying the influence of particle breakage on the critical state behaviours of granular soils. For instance, Bandini and Coop [4] performed triaxial tests on carbonate Dog's Bay sand and found that the change of grading due to particle breakage causes both the vertical and rotational movements of the CSL. Yu [92] conducted a series of triaxial tests on natural and pre-crushed sands. The author revealed that the CSL evolved in a combination of a downward translation and an anticlockwise rotation in the e -logp 0 plane as particle crushing increases, while in the p 0 -q plane, a two-segment linear CSL can be captured. Zhang and Luo [97] carried out triaxial tests on a dense calcareous sand. They found that in the p 0 -q plane, the CSL exhibited a downward shift due to particle breakage, resulting in a change of the critical state void ratio of the material.
The breakage of rockfill materials at high stress levels has also been a concern for geotechnical projects. Xiao et al. [86] investigated the influence of particle breakage on the position of CSL of Tacheng rockfill material using large-scale triaxial tests. The results show that the critical state stress ratio in the p 0 -q plane is not affected by particle breakage, but the CSL in the e -logp 0 plane shifted downward and rotated counterclockwise as particles break. Based on the experimental observations, Xiao and Liu [85] proposed an elastoplastic constitutive model for rockfill materials considering particle breakage, which can well predict the stress-strain and particle-breakage behaviours of rockfill materials at various confining pressures. Guo et al. [23] investigated the effect of particle breakage on the critical state strength of rockfill materials using large-scale triaxial tests and found a passive correlation between critical state stress ratio and particle breakage ratio, rather than being a constant. This is different to the findings by Xiao et al. [86].
However, in laboratory tests, particle breakage is generally characterized using pre-crushed sands with different initial gradings [4,77,79,92], or comparing the soil's preshearing and post-shearing gradings [23,83,86]. The main shortcoming of these approaches is that the continuous evolutions of soil's grading and fabrics during shear cannot be studied.
Discrete element modelling (DEM) is an effective way to capture the changes in microscopic fabrics within the soil samples during shear [19,38,54,56,64,96]. It provides a good channel to understand the micromechanical behaviours of granular materials. As to the simulations of particle breakage, the information involving the evolutions of particle size and fabric anisotropy can be easily obtained using DEM [37,44,70,101]. Hanley et al. [25] performed triaxial compression simulations on the assemblies of crushable spherical particles. The results show that in the p 0 -q plane, a curved envelope of peak shear resistance can be attained when particle breakage is considered, but the CSLs are both linear whether particle breakage is considered or not. Ciantia et al. [13] simulated triaxial tests using DEM to study the effect of particle fragmentation on CSL. The results demonstrate that the position of CSL is strongly related to the breakage-induced grading evolution. Shi et al. [69] studied the coupled effects of intermediate principal stress ratio b and particle breakage on the critical state behaviours using a series of true triaxial DEM simulations on Nansha calcareous sand. The authors found that particle breakage increases with the increase of b, leading to the decrease of both the intercept and the slope of the CSL in the e-p 0 plane.
Existing studies have been focused on the effects of particle breakage on the critical state behaviours of granular soils under drained conditions but little under undrained conditions. Furthermore, the findings from the limited literature regarding the effect of drainage conditions on the CSLs are still not concordant. Some researchers [5,6,24,55] stated that the CSL is independent of the drainage condition, but others found a difference in the positions of CSLs in drained and undrained triaxial tests [88]. However, all the above discussions are in the absence of particle breakage. Heitor et al. [28] performed drained and undrained triaxial tests on coal wash and found that the incidences of particle breakage under drained and undrained shear conditions are different, which in turn strongly affects the critical state of the material. Yu [93] examined the effects of particle breakage on the critical state of sands under undrained shear conditions using triaxial tests. The authors found that under the same initial void ratio, the p 0 values at the critical state decrease as particles break in the e -logp 0 plane, whilst in the p 0q plane, the critical state of the soil moves downward along the CSL as particle breakage increases.
In this study, triaxial tests under undrained shear conditions are simulated in DEM using the constant volume method in 3-dimensional particle flow code (PFC3D). The model of single particle crushing is established via the octahedral shear stress criterion and the crushing mode conforming to an Apollonian sphere packing. The parameters for particle breakage are obtained from published experimental results on Nansha calcareous sand. The effects of particle breakage on the undrained shear behaviours of granular soils with different initial void ratios and confining pressures are investigated, with emphasis on the variation of CSLs in both p 0 -q and e-p 0 spaces. The relationship between the macro-and microbehaviours of the soil assemblies is also examined by studying the evolutions of the coordination number, redundancy ratio and overall anisotropy of soil fabrics during the undrained shear.

DEM modelling
2.1 Single-particle crushing model When DEM is applied to simulate the particle breakage, two common solutions have often been used. One is the bonded particle method [7,39,47,63], in which porous agglomerates are created and can later disaggregate when the contact force (stress) acting on the bonding point exceeds the assigned bonding strength. Another is the particle replacement method [2,33,45,75,101], in which once the tensile stress within a particle induced by the contact force exceeds the breakage strength, the original particle breaks and is replaced with an equivalent group of subparticles. It is worth noting that both methods rely upon the use of spheres. The main weakness of the bonded particle method is the low computational efficiency due to the large number of elementary particles required for the generation of agglomerates. Therefore, the particle replacement method is adopted here.
When using the particle replacement method to establish the single-particle crushing model, particle breakage criterion and fracture mode are two key issues. The soundness of different breakage criteria has been discussed by de Bono and McDowell [16]. They recommended that the octahedral shear stress criterion is an optimal choice to consider the breakage of sand particles under multiple contacts, which can be expressed as: where r 1 , r 2 and r 3 are the major, intermediate and minor principal stresses acting on the particle, respectively, which can be calculated from the stress tensor: where V is the bulk volume; N c is the total number of contacts of the particle; f j c is the jth component of the contact force at the contact point; l i c is the ith component of the branch vector of the contact.
According to McDowell and de Bono [50], when a particle under multiple contacts is considered, the octahedral shear strength q n can be calculated as: where d is the particle diameter; d 0 is the characteristic particle diameter; r 0 is the characteristic tensile stress at which 37% of the particles survive; m is the Weibull modulus. The detailed definitions regarding characteristic tensile stress and survival probability can be found in Weibull [81] and Nakata et al. [51]. In this model, particle breakage occurs when the induced octahedral shear stress q oct is greater than or equal to the octahedral shear strength q n of the particle (i.e. q oct C q n ).
As to the selection of fracture mode, the number of spawned fragments and the mass conservation are two important problems that need to be concerned. As suggested by Zhou et al. [101], the present fragment replacement modes can be divided into three categories: major splitting, substantial crushing and the Apollonian sphere packing modes. Due to the particle breakage and fragmentation, the particle size distribution of granular soils is expected to have an ultimate fractal distribution of 2.5 * 2.7 [1,18,49,76]. The Apollonian sphere packing is a typical fractal packing mode with a fractal dimension of 2.47 [8], which is quite suitable for imitating the fractal distribution of fragments. Ciantia et al. [12] investigated the computational efficiency and accuracy using different numbers of fragments satisfying the Apollonian sphere packing, in which the fragment numbers of 2, 3, 4, 5, 8, 14 and 57 are examined, and the authors recommended that an Apollonian packing of 14 fragments is an ideal choice to balance the computational cost and accuracy in DEM modelling. Therefore, in this study, a fracture mode similar to that introduced by Ciantia et al. [12] is adopted. As shown in Fig. 1a, a clumped particle consisting of 14 elementary sub-spheres with the distribution conforming to an Apollonian packing is generated first in PFC3D. The particle size of the clumped particle is decided by the diameter of the circumscribing sphere. When the induced octahedral shear stress acting on the clumped particle satisfies the predefined particle crushing criterion as described in Eq. (3), the clumped particle separates into 14 spherical fragments. Meanwhile, multigenerational crushing is permitted as shown in Fig. 1b. The spherical particle generated in the first generation can be further split into 14 new small fragments according to the same fracture mode. It should be indicated that as to the breakage of the original clumped particles, the mass conservation is maintained. The mass loss only happens when the fracture of the child particles is considered. Although the 47% mass loss occurs when a child particle breaks into a group of 14 fragments obeying the Apollonian packing [12], the total mass loss of the sample is negligible (less than 2%) compared with the total mass of the original particle system, due to the low proportion of the failed child particles in the simulations. This will be detailed later.

Sample preparation
To prepare the samples, a polydisperse assembly of 1906 clumped particles (i.e. 26,684 elementary sub-spheres) with random orientations was generated in PFC3D. The particle size of the assembly ranges from 0.25 to 2 mm, which is the same as that of Nansha calcareous sand [9]. The particle size distributions of the DEM samples and Nansha calcareous sand are compared in Fig. 2. Six rigid and frictionless walls are built to form a cubic triaxial box with an internal dimension of 5.5 mm. The particles are placed into the box, and then, the six walls move inward at 0.05 mm/s until the target initial void ratio e 0 is reached. Three initial void ratios, i.e. e 0 = 1.0, 1.1, 1.3, are considered to study the effect of initial density on the behaviour of the soil. It is worth noting that the initial void ratio in this analysis is calculated using the solid volume and the total volume of the numerical sample to avoid the potential fluctuation of the void ratio, resulting from the release of the interior voids due to the breakage of the clumped particles. The same strategy has also been used by Xu et al. [84] in their DEM study. The main shortcoming of this definition is that the calculated void ratio is a little greater than that of physical specimens with the same particle size distribution. After the sample reaches the target e 0 , isotropic confining stress r 3c is applied on the sample by further moving the boundary walls inward. Five confining stress levels, i.e. r 3c = 50, 100, 200, 300, 400 kPa, are used in the simulations. It should be indicated that small volume shrinkage occurs when the sample is subjected to confining stress. A representative DEM sample after isotropic compression is shown in Fig. 3a.
A linear contact stiffness model with Coulomb sliding criterion was used to describe the inter-particle contact behaviour. In addition, considering the distinct angular shape of calcareous sand particles [77,79], a linear rolling resistance model was introduced to overcome the drawback of approximately spherical clumped particles used in the model. The local non-viscous damping proposed by Potyondy and Cundall [61] was utilized to ensure the quasi-static equilibrium of the particle system and to consider the existence of non-rattler particles, i.e. particles with zero contacts [36]. The micromechanical parameters used in the model are summarized in Table 1. The parameters for particle breakage criterion as described in Eq. (3) are from the single-particle crushing test results on Nansha calcareous sand reported by Ma et al. [46].
The model verification and parameter calibration were accomplished by comparing the numerical modelling with the triaxial drained shear test results on Nansha calcareous sand [9]. The initial void ratios of the sand sample and DEM sample are 0.931 and 0.981, respectively. The confining stresses for the two samples are 100 kPa. Figure 4 shows the comparison of the shear strength and volume deformation responses between the calcareous sand and DEM samples. The influence of the rolling resistance coefficient l r on the macro-behaviours of DEM samples is also demonstrated in the figure. As shown in the figure, the strength at 25% of strain and the volume deformation response of the DEM sample (l r = 0.5) are both in line with those of the calcareous sand, though the peak strength (a)

(b)
Break into Fig. 1 The particle crushing mode: a the original clumped particle; b the fracture mode of the child fragments Fig. 2 Comparison of particle gradings between the DEM samples and the Nansha calcareous sand [9] of the DEM sample (l r = 0.5) is a little lower than that of the calcareous sand. It can also be seen that when the contact stiffnesses and the inter-particle friction coefficient are kept constant, the peak shear resistance and volume dilation both increase with l r . The optimum choice for l r is 0.5 for the present calibration.

Simulation of undrained triaxial shear tests
To simulate an undrained shear condition in DEM, two approaches are commonly used. One is DEM coupled with computational fluid dynamics (CFD) [11,42,59], in which the kinematics of the porous fluid is solved using CFD and the movement of the particles is calculated by DEM. The advantage of the CFD-DEM coupled method is that the interaction between fluid and particle can be analysed at the particle level, yet the main shortcoming is the low computational efficiency. Another approach is the constant volume method [21,34,67,71], in which the volume of a sample is kept constant during shear to simulate an undrained condition. Although the interaction between fluid and particles cannot be directly studied using the constant volume method, the excess pore water pressure due to compression or dilation of the samples during shear  can be easily obtained. The validity of the constant volume method has already been verified by Shafipour and Soroush [67], and Zhang and Rothenburg [98]. They stated that compared with the CFD-DEM method, the constant volume method can capture the same undrained shear responses of granular soils with less computational time. In this regard, the constant volume method is utilized here. After the application of the isotropic confining pressure, the servo-control mechanism is turned off in PFC3D. To achieve a constant volume under the undrained triaxial shear condition, during the shear stage, the top and bottom walls move inward at a constant velocity with v z-= 0.05 mm/s, and the lateral walls in both x and y directions move outward at the velocity of: where l x and l z are the real-time sample lengths in the x and z directions during shear as shown in Fig. 3b. According to Eq. (4), the initial values of v x and v y are 0.025 mm/s. The lateral velocities update every 50 timesteps to assure the constant volume of the sample during shearing (the maximum variation of sample volume is within 0.002%). The shear stage terminates when the major principal strain e 1 (i.e. axial strain in the z direction) reaches 40%, where the strain is always calculated using the initial size of the samples. In the simulations, grain crushing is only considered during the shear stage. As a comparison, simulations were also conducted on uncrushable particles under the same conditions to examine the effect of particle breakage on the critical state behaviours of the sample.

Strength and excess pore water pressure characteristics
To discuss the micromechanical behaviours of the DEM samples under undrained shear conditions, the mean effective principal stress p 0 , generalized shear stress q and excess pore water pressure u defined in the following equations are used: where r 0 1 , r 0 2 and r 0 3 are major, intermediate and minor effective principal stresses, respectively. The effective stress, taking r 0 1 for example, is calculated as the average value of the effective stresses acting on the top and bottom walls in the z direction, as shown in Fig. 3. The effective stress on a wall is calculated as the total normal contact force acting on the wall divided by the real-time wall area. As depicted in Eq. (7), the excess pore water pressure can be estimated as the difference between the initial confining stress and the real-time effective stress acting on the lateral walls. Figure 5a-c shows the variation of the generalized shear stress ratio (q/p 0 ) with the major principal strain (e 1 ) of the DEM samples with different initial void ratios (i.e. e 0 = 1.0, 1.1 and 1.3). To effectively investigate the effect of particle breakage on the shear strength of the sample, the  Fig. 5 that the peak value of q/p 0 , i.e. (q/p 0 ) p , declines due to particle breakage during shear, and the effect of particle breakage on (q/p 0 ) p is more significant for denser soils. Taking r 3c = 50 kPa for example, when e 0 = 1.3 (Fig. 5c), the value of (q/p 0 ) p declines from 1.54 to 1.51 due to particle breakage, and the difference is 0.03. As e 0 decreases to 1.0 (Fig. 5a), the value of (q/p 0 ) p reduces from 1.66 for the uncrushable sample to 1.54 for the crushable one, with a drop of 0.12 due to particle breakage. The reason for this discrepancy is caused by the different extents of particle breakage in the samples with different initial densities. This will be further discussed in the following sections.
It is hard to obtain reliable results at strain levels greater than 40% due to the limitation of the strain levels that can be simulated in DEM. As shown in Fig. 5, some small fluctuations still exist in the stress ratios at the strain level of 40%. Therefore, at this stain level the samples may not be strictly at the critical state. Considering that many researchers have adopted the 40% strain level as the critical state [53,100], the state of the samples achieved in the modelling is assumed at an approximate critical state. In the following discussions, the term critical state is used to describe the behaviour of the samples at the strain level of e 1 = 40%. It is found in Fig. 5 that particle breakage has little impact on the critical state strength of the samples under the undrained shear condition, irrespective of initial void ratio or confining pressure. Compared with the influence of e 0 , the effect of r 3c on the undrained shear strength of the samples is less, which is independent of particle breakage. Figure 6a-c depicts the variation of excess pore water pressure (u) with the major principal strain (e 1 ) of DEM samples with e 0 = 1.0, 1.1 and 1.3, respectively. In Fig. 6, the positive values of u indicate the volume compression tendency of the samples under the undrained shear, and the negative values correspond to the volume dilation tendency of the samples. It can be seen in Fig. 6 that although the variation of u with strain follows a similar trend for crushable and uncrushable samples, greater positive values and smaller negative values of u can be captured in soils with crushable particles, irrespective of e 0 and r 3c . This indicates the greater contraction tendency of the samples due to particle breakage, which agrees with the experimental findings in the undrained shear tests by Yu [94]. The effect of e 0 on the variation of u can also be found in Fig. 6: the higher the e 0 , the higher the positive value of u. As shown in Fig. 6c, for the crushable samples with e 0 = 1.3, the maximum value of u almost reaches the initial confining pressure for all the confining stress levels. This means that the sample reaches the state of initial liquefaction. Figure 6 also demonstrates the effect of r 3c on the variation of u. It can be seen in the figure that r 3c only influences the value of u but not the variation tendency of u with e 1 . Higher positive and lower negative values of u are exhibited with higher r 3c , especially for the sample with crushable particles.  Fig. 7 that due to the increment of u in the initial shear stage as shown in Fig. 6, the mean effective principal stress p 0 first decreases from the initial confining pressure to the phase transformation point (PTP) for both crushable and uncrushable samples. Compared to that of the uncrushable assemblies, there are sudden drops of p 0 values in the crushable assemblies at the start of the shearing. This is due to the increase of u (shown in Fig. 6) induced by particle breakage as discussed later, and the higher the confining pressure, the greater the drop in p 0 . The figure also shows that after PTP, the value of p 0 gradually increases along an approximate oblique line (i.e. shear failure envelope) for both crushable and uncrushable samples. But the slopes of the failure envelope, i.e. M f , are different for crushable and uncrushable samples, especially for the denser samples. Taking e 0 = 1.0 for example, the M f decreases from 1.62 to 1.50 when particle breakage is considered. According to Wautier et al. [80] and Yilmaz et al. [91], the fundamental reason for the decrease in M f with particle breakage can be attributed to the addition of fine particles within the particle assembly, which suppresses the volume dilation and reduces the shear resistance of the sample.

Evolution of particle breakage
The relative breakage ratio B r proposed by Hardin [26] is used to quantitatively study the evolution of particle breakage during shear: where B p is the initial breakage potential, which can be expressed as the area surrounded by the initial grading curve and the line of 0.074 mm grain size, as shown in Fig. 8; B t is the breakage potential represented as the area enclosed inside the current grading, initial grading and the 0.074 mm grain size line. Figure 9 shows the variation of B r of the three samples during shearing. The figure shows that, for all the samples, particle breakage initiates as soon as shearing starts. At lower strain levels, i.e. less than 2.5%, the rate of particle breakage increases as the soil becomes denser, and the higher the confining pressure, the higher the rate of particle breakage. For the two samples with e 0 = 1.0 and 1.1, B r increases sharply at the start of the shearing under higher confining pressure levels, e.g. r 3c C 300 kPa. Nearly 50% of the total particle breakage finished within 2.5% of strain. This results in the sudden increase of pore water pressure at the initial shear stage of the samples as shown in Fig. 6. As strain increases, the rate of particle breakage decreases. At larger strains, i.e. [ 20%, the relative breakage ratio increases almost linearly at a much lower rate. The figure also shows that the relative breakage ratio at the critical state (e 1 = 40%), (B r ) cs , is much affected by confining pressure for the looser sample (e 0 = 1.3), but is less sensitive to confining pressure for the two denser samples, which can be better seen in Fig. 10. This suggests that fewer particles will break in loose samples under lower confining pressures.
The evolution of particle breakage can also be described by the crushing probability, P c , which can be calculated as the ratio of the crushed particle number to the total number of the specific generation. A typical case of a dense sample (e 0 = 1.0) under high confining stress (r 3c = 400 kPa) is selected to show the evolution of P c under the undrained shear loading. As depicted in Fig. 11, the whole crushing process is mainly dominated by the fracture of the original clumped particles, whilst the crushing probability of the child fragments is very low. When the strain reaches 40%, the P c of the original clumped particles is 9.3%, i.e. 177 original particles crushed among the 1906 total original particles. As stated earlier, no mass loss occurs in this crushing stage. The P c of the child fragments is 0.65% at e 1 = 40%, i.e. 16 child particles crushed among the 2478 total generated child particles, which leads to 1.8% mass loss of the whole particle system. No breakage of the grandchild particles has been observed in the simulations due to the high breakage strength of fine particles.

Critical state behaviours
In this section, the effects of particle breakage and drainage condition on the CSLs in both p 0 -q and e-p 0 spaces are discussed. To investigate the effect of a drainage condition, drained triaxial tests on the samples with the same e 0 and r 3c conditions are also simulated using both the crushable and uncrushable particles. Figure 12 illustrates the p 0 cs -q cs (the p 0 and q at the critical state, respectively) relationship of the samples at the critical state. The abbreviations of CU and CD in the figure depict the undrained and drained triaxial simulations, respectively. It can be seen in Fig. 12 that, at the critical state, all the samples follow nearly the same Mohr-Coulomb failure criteria regardless of drainage condition or crushability of the particles: where M cs is the critical state stress ratio, which is about 1.38. The value is close to the laboratory test result (1.44) reported by Yu [92]. This confirms that the critical state Fig. 8 The definition of relative breakage ratio B r [26] stress ratio of granular soils is independent of the initial state, drainage condition and particle breakage.
As suggested by Li and Wang [40], the CSLs are plotted in the e cs -(p 0 cs /p a ) n plane and fitted with the following equation: where e cs is the void ratio at the critical state, e C and k are the intercept and the slope of the CSL, respectively; p a is the atmospheric pressure (taken as 101 kPa here); n is the material coefficient and often taken as 0.7 for sand as suggested by Li and Wang [40], and Hanley et al. [25]. Figure 13 represents the e cs -(p 0 cs /p a ) n relationship. It is to be noted that, in the figure, the void ratios at the critical state are slightly different to the initial void ratios. This is due to the different initial deformations caused in the assemblies under different confining pressures. So even for the samples with the same initial void ratio, the void ratios at the critical state are different under undrained shear conditions. It can be seen that the crushable and uncrushable samples follow different relationships. The critical state line moves downward and rotates clockwise when particle breakage is considered. The comparison shows that the drainage condition has little effect on the CSL of the samples.
For crushable samples, Fig. 14a illustrates the relationship of (B r ) cse csp 0 cs . It can be seen in the figure that a unique correlation exists between (B r ) cs , e cs and p 0 cs , regardless of the initial state and drainage condition. A similar relationship has been observed on rockfill materials considering particle breakage reported by Xiao et al. [86], and Xiao and Liu [85]. Figure 14b plots the relationship between (B r ) cs and (p 0 cs /p a ) 0.7 of all the samples. It can be seen that (B r ) cs increases nearly exponentially with (p 0 cs / p a ) 0.7 . Meantime, the increasing tendency of (B r ) cs exhibits two-stage characteristics, when (p 0 cs /p a ) 0.7 C 4, the increasing ratio of (B r ) cs declines, indicating that particle breakage gradually enters into a diminishing stage with increasing stress level. Figure 14c shows the projection of the data points to the e cs -(B r ) cs plane. It can be seen that although e cs decreases with the increase of (B r ) cs , the relationship is not linear, and an exponential fitting curve can be used to express the e cs -(B r ) cs relationship.

Coordination number
The variation of contact number between particles has been used to describe the evolution of soil volume. The average contact number within a particle system can be represented with the coordination number Z: where N p is the total number of particles, which needs to be updated during the calculations; N c is the total number of contacts of the particle system. Figure 15a-c depicts the variation of Z with the major principal strain e 1 of the samples under undrained shear. It can be seen that under the undrained shear condition, the Fig. 10 The relationship between (B r ) cs and r 3c Fig. 11 The evolution of crushing probability (e 0 = 1.0, r 3c = 400 kPa) Fig. 12 The CSLs in the p 0 -q plane Fig. 13 The CSLs in the e-p 0 plane variation of Z with e 1 agrees well with that of u shown in Fig. 6. For all the samples, the samples tend to contract at the start of the shearing. So the soil particles tend to move into unoccupied voids, leading to the drop of the contact numbers and the drop of effective stress, which in turn results in the rise of pore water pressure (Fig. 6). In this explanation, the rise of pore water pressure is explained as the consequence of the drop of effective stress rather than the other way around as observed in laboratory tests. This is due to the constant volume shearing technique used in the modelling of undrained shear. At larger strain levels, i.e. [ 20%, the average contact number in the uncrushable samples keeps increasing, but stabilizes in the crushable samples. This explains the increase of negative pore water pressure in the uncrushable samples, indicating an expanding tendency of the samples. In the crushable samples, due to the increase of B r with strain as shown in Fig. 9, the average contact number remains almost constant at strain levels greater than 20%. This greatly reduced the increasing tendency of average contact numbers. It is worth noting that the average contact number of the crushable samples with e 0 = 1.3 drops almost to zero at the start of

Redundancy ratio
As suggested by Kruyt [36], and Pouragha and Wan [62], the redundancy ratio R f is an important indicator to describe the inter-particle contact behaviour and the stability of a particle system under shear loading. The redundancy ratio can also be investigated by comparing the number of constraints at contacts with the degrees of freedom of the particles. When the Coulomb sliding criterion and the rolling resistance model are introduced to describe the contact behaviour between two particles, the concepts of sliding contact, rolling contact and slidingrolling contact are used to describe the failure modes of the contact between the two particles: sliding failure, rolling failure and sliding and rolling failure, respectively. As proposed by Zhou et al. [100], the redundancy ratio R f can be expressed as: where R nsr is the ratio of the number of no sliding nor rolling contacts to the total number of contacts; R s , R r and R sr are the ratios of the number of sliding contacts, rolling contacts and sliding-rolling contacts to the total number of contacts, respectively. According to Zhou et al. [100], for an equilibrium state or a solid-like particle system, R f should be greater than 1, and R f = 1 is the boundary between solid-like and liquid-like behaviour. Figure 16a-c shows the variation of R f with e 1 of the samples under the CU shear test. It can be seen that the variation of R f with strain is similar to that of the coordination number Z. This indicates that the coordination number plays a dominant role in the evolution of contact behaviour. Compared with uncrushable samples with the same e 0 and r 3c , the R f significantly decreases when particle breakage is considered. For crushable samples with large e 0 , e.g. e 0 = 1.3, the value of R f is even less than 1.0 in the initial shear stage, which is resulted from the initial liquefaction of the samples under the undrained shear.
The critical state redundancy ratio, i.e. (R f ) cs , is extracted from Fig. 16 at e 1 = 40% to analyse the relationship between the microstructure and macroscopic behaviours of the sample. Figure 17a shows the relationship of (R f ) cse csp 0 cs in a 3D space. The results from the CD simulations are also plotted in the figure as a comparison. It can be seen in Fig. 17a that for both crushable and uncrushable samples, the data points exhibit a good correlation, which is independent of drainage conditions. To reveal this correlation more clearly, 2D projections to the (R f ) csp 0 cs and (R f ) cse cs planes are shown in Fig. 17b and c, respectively. As shown in Fig. 17b, the relationship between (R f ) cs and (p 0 cs /p a ) 0.7 is linear for both crushable and uncrushable samples, irrespective of initial density and drainage condition. For crushable samples, (R f ) cs increases with (p 0 cs /p a ) 0.7 at a much lower rate, i.e. about half the rate of that of uncrushable samples. For the samples with the same e 0 , the values of (R f ) cs and (p 0 cs /p a ) 0.7 are much smaller (about half or less) than those of uncrushable samples. The main reason is that compared with the uncrushable sample, particle crushing significantly decreases the average contact number within the particle system, as shown in Fig. 15. Similar to the (R f ) csp 0 cs relationship, a linear regression can also be obtained to express the relationship between R f and e cs for both crushable and uncrushable samples. It is well known that at the macro-level, the critical state of the soil can be explained as a constant volume state under constant stress. The present DEM findings on the unique (R f ) cse cs p 0 cs relationship prove that at the micro-level, the critical state corresponds to a specific static equilibrium of internal structure of the soil, regardless of particle breakage or drainage condition, and this static equilibrium can be quantitatively characterized by the critical state redundancy ratio. The findings provide solid evidence for the establishment of a micromechanical constitutive model for describing the critical state of the soil.

Fabric and force anisotropies
In addition to the coordination number and redundancy ratio, the evolution of fabric and force anisotropies due to shear loading is another important micromechanical phenomenon that needs to be concerned. Earlier micromechanical studies verified that the soil fabric at the critical state exhibits a clear anisotropic structure [22,48,58]. The degree of fabric and force anisotropies can be quantified by the following anisotropy coefficients [69,103]: ; a n ¼ ffiffiffiffiffiffiffiffiffiffiffiffi ffi 3 2 a n ij a n ij r ; and a t ¼ where a c , a n and a t are the anisotropy coefficients of contact normal, normal contact force and tangential contact force, respectively; a c ij , a n ij and a t ij are the second-order anisotropy tensors of contact normal, normal contact force and tangential contact force, respectively. Figure 18a-c shows the evolution of the overall anisotropy (a c ? a n ? a t ) of the samples. It is found that the evolution of overall anisotropy is similar to the variation of stress ratio shown in Fig. 5. As depicted in Fig. 18, significant soil anisotropy develops with strain. Taking the dense sample (e 0 = 1.0) with crushable particles for example, the peak value of the overall anisotropy reaches the range of 3.85 * 3.95. At the critical state, the crushability of the samples and confining pressure have little impact on the overall anisotropy of the samples. This can be better seen in Fig. 19, which shows the relationship between the critical state overall anisotropy (a c ? a n-? a t ) cs and the critical state stress ratio (q/p 0 ) cs of the samples. The figure shows a linear relationship between the critical stress ratio and the overall anisotropy. This suggests that at the critical state, the ultimate soil strength is determined by the extent of fabric and force anisotropies, and the overall anisotropy ratio could be directly related to the critical state friction angle, regardless of particle breakage and drainage condition.

Conclusions
A number of DEM simulations were conducted on the assemblies of crushable and uncrushable particles under drained and undrained triaxial shear conditions, to study the effects of drainage conditions and particle breakage on the critical state behaviours of granular soils. In the analysis, the samples are assumed to reach their critical states at the axial strain level of 40%. The variation of shear resistance, excess pore water pressure, the interaction between particles, and the evolution of soil fabrics are studied to reveal the relationship between microscale parameters and macroscale soil behaviour. Based on the limited number of DEM simulations, it is found that under the undrained condition: (1) the crushing of particles could reduce the peak shear stress ratio, especially for denser samples. The confining stress level has little impact on the peak shear stress ratio. The crushability of the particles would not greatly affect the critical shear stress ratio.
(2) Crushing of particles would increase the excess pore water pressure accumulated at the initial shear stage and decrease the negative pore water pressure in the subsequent loading stage. This is due to the breakage of the particles, which increases the shrinkage and decreases the expansion tendency of the assemblies. For looser samples, the excess pore water pressure could be as high as the confining pressure, which indicates the initial liquefaction of the samples. The relative breakage ratio (B r ) is affected by the initial void ratio (e 0 ), and the smaller the e 0 , the greater the B r . This explains the greater impact of particle breakage on the critical state stress ratio of denser samples. The relative breakage ratio at the critical state, (B r ) cs , has a near positive linear relationship with confining pressure.
It is also found that: (1) the particle breakage and drainage condition do not affect the slope of the critical state p 0 -q line. This again confirms that the critical state of a soil is an intrinsic property of the soil and is not affected by the initial state, drainage condition and crushability of the particles. (2) In the e cs -(p 0 cs /p a ) 0.7 space, the drainage condition does not affect the critical state line. The critical state line translates downward and rotates clockwise when particle breakage is considered.
(3) At the critical state, the relative breakage ratio (B r ) cs , the mean effective stress (p 0 cs ), and void ratio (e cs ) are correlated. There is an exponential relationship between (B r ) cs and (p 0 cs /p a ) 0.7 , and e cs and (B r ) cs . (4) The redundancy ratio R f reduces as particles break.
At the critical state, there is a unique relationship between (R f ) cs , e cs and (p 0 cs /p a ) 0.7 , regardless of the crushability of the samples and the drainage condition. This indicates that the critical state redundancy ratio can be an indicator of the internal stability of soil structures. (5) The critical state overall anisotropy (a c ? a n ? a t ) cs is an important indicator of the shear strength parameter. At the critical state, there is a unique linear relationship between (a c ? a n ? a t ) cs and the critical state stress ratio (q/p 0 ) cs of the samples, regardless of the crushability of soil particles and the drainage condition. Fig. 19 The relationship between the overall anisotropy and the generalized shear stress ratio at the critical state Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.
Funding Open Access funding enabled and organized by CAUL and its Member Institutions.