Micromechanical inspection of incremental behaviour of crushable soils

In granular soils grain crushing reduces dilatancy and stress obliquity enhances crushability. These are well-supported specimen-scale experimental observations. In principle, those observations should reflect some peculiar micromechanism associated with crushing, but which is it? To answer that question the nature of crushing-induced particle-scale interactions is here investigated using an efficient DEM model of crushable soil. Microstructural measures such as the mechanical coordination number and fabric are examined while performing systematic stress probing on the triaxial plane. Numerical techniques such as parallel and the newly introduced sequential probing enable clear separation of the micromechanical mechanisms associated with crushing. Particle crushing is shown to reduce fabric anisotropy during incremental loading and to slow fabric change during continuous shearing. On the other hand, increased fabric anisotropy does take more particles closer to breakage. Shear-enhanced breakage appears then to be a natural consequence of shear-enhanced fabric anisotropy. The particle crushing model employed here makes crushing dependent only on particle and contact properties, without any pre-established influence of particle connectivity. That influence does not emerge, and it is shown how particle connectivity, per se, is not a good indicator of crushing likelihood.

Two well-established observations arising from earlier research are that (a) soil crushability is enhanced by stress obliquity, and (b) plastic deformation in crushable soils is associated with less dilatancy, or, equivalently, that plastic flow is more volumetric when crushing is present. These specimenscale observations have been reproduced by several DEM models of crushable soils [13,15,34]. Specimen-scale responses of soils may be considered as emerging responses, i.e. responses that reflect underlying micromechanical interactions at the grain or element scale. This general principle has not been yet systematically applied to crushable soils, and several aspects of their micromechanics remain relatively obscure.
Laboratory micromechanical investigations of breakable soils are possible [1] but very difficult. On the other hand, DEM simulation offers readily accessible information at the element scale, which may be used to uncover relevant micromechanics. Micromechanical analysis of breakable DEM models has mostly focused on the conditions that lead to particle breakage and the nature of the obtained fragments [27,29,38,55]. Less attention has been paid to the underlying mechanisms that, for instance, may explain the more volumetric nature of crushing-induced plastic flow. This bias is partly due to the use of crushable agglomerates to represent breakable particles in DEM models. When using agglomerates, it is difficult to track inter-granular fabric evolution [19]. Agglomerates may be too compliant [6] or have multiple contact points [64], making fabric measures difficult to interpret. Modelling particles as agglomerates also carries significant computational costs; this limits the number of particles in simulation and, therefore, the robustness of micromechanical statistics. Still, a general observation derived from crushable agglomerate studies [6,61,64] is that the development of contact force anisotropy during shear is strongly affected by crushing.
To avoid the problems associated with agglomerates, crushable DEM models based on the splitting particle concept [41] may be used. Zhou et al. [68] used a splitting particle model to study fabric evolution during true-triaxial simulation of breakable particles. However, the study focused on the effects of intermediate principal stress ratio and the connection between breakage and fabric evolution was not investigated.
Ciantia et al. [16] used a DEM model based on the splitting particle concept to explore the incremental response of crushable soils. They showed that grain crushing causes the plastic flow to change direction and that this direction is unique and independent on the loading direction. Crushinginduced plastic flow was observed to be more volumetric than when crushing is inhibited. However, the micromechanics of that response was not investigated.
This study aims to fill that gap, by investigating how particle breakage modifies the micromechanics underlying the incremental stress response of granular specimens. Response envelopes were originally proposed by Gudehus [32] as a means to classify constitutive models. Incremental response envelopes [24,25,54] reveal important characteristics soil behaviour such as incremental nonlinearity, associative behaviour, uniqueness of flow rule, etc. Measuring incremental responses in the laboratory is difficult, [21,22], because of the highly accurate strain measurements required and the need to test several identical samples (one for each probe). For this reason DEM models have been frequently used to study the incremental behaviour of granular materials [3,7,11,47,58,60].
Taking advantage of the efficiency of the DEM crushing model of Ciantia et al. [16], results have been obtained from a large number of simulations ([ 350). In what follows the DEM model is first described; the methods applied for micromechanical investigation are then introduced before presenting some relevant numerical results. The focus of the paper is to examine the relationship between crushing and the features and evolution of the microstructure. It also aims to explain mesoscale observations such as the enhanced volumetric component of crushing-induced strains, the uniqueness of plastic flow direction and the cascading nature of particle breakage.

DEM particle crushing model
The DEM model for crushable soil employed here was described by Ciantia et al. [15]. It was implemented as an userdefined model in the PFC3D code [35] and for completeness an overview is given here. Coulomb friction and the simplified Hertz-Mindlin contact model are used to describe the interaction of spherical elements. Particle rotation was inhibited to give a good match to experimental data without using nonspherical particles. This approach, successfully applied in previous research [2,8,9,59], is equivalent to applying a very large value of rotational resistance at the contacts.
The failure criterion is based on work by Russell and Muir Wood [48] and Russell et al. [49]; a two-parameter material strength criterion is used together with consideration of the elastic stresses induced by point loads on a sphere. A particle subject to a set of external point forces reaches failure when the maximum applied force reaches the following limit condition: where r lim is the limit strength of the particle and A F is the contact area. To incorporate the natural material variability into the model, the particle limit strength, r lim , is assumed to be normally distributed for a given sphere size. The coefficient of variation of that distribution, var, is taken to be a material parameter. Following McDowell and De Bono [43], the mean strength value for given sphere size, r lim , depends on the particle diameter: where r lim 0 is a material constant and f(var) indicates the effect of variability of particle strenght. A F depends on the contact force and the particle's elastic properties; applying Hertzian contact theory, the following expression for the breakage criterion is obtained: where r 1 and r 2 are the radii of the contacting spheres and E i , m i are the Young's moduli and Poisson's ratio, respectively. Note that this breakage criterion does not involve exclusively the maximum force on the particle: there is a strong inbuilt dependency on the characteristics of the contacting particles.
Once the limit condition is reached, the spherical particle will split into smaller inscribed tangent spheres. The crushed fragments assume the velocity and material parameters of the original particle apart from the intrinsic strength (r lim0 ) which is randomly assigned sampling its normal distribution. Ciantia et al. [15,16] concluded that a 14-ball crushed configuration can adequately represent macroscopic behaviour. The capability of this model to capture real test behaviour has been demonstrated by Ciantia et al. [14,17].
Ciantia et al. [16] calibrated this model to reproduce the observed experimental behaviour of Fontainebleau sand at low pressures ( Fig. 1). Particle-scale information, from single particle crushing tests on quartz sands, was also used to select crushing model parameters. The grain size distribution is also close to that of Fontainebleau sand [66]. More recently, Ciantia et al. [17] refined this calibration, using new data to obtain better experimental agreement at very high pressures. The work presented here is not aimed to match laboratory experiments and, for coherence with previous work on response envelopes, the model parameters applied ( Table 1) are those of Ciantia et al. [16].

Numerical specimen
A DEM specimen was created by filling a cube with spheres with particle diameters ranging from 0.1 to 0.4 mm. The cube side lengths of 4 mm (Fig. 1a) contained close to 10,000 particles. This is a representative volumetric element (REV), as small as possible to ease computation, but big enough so that the observed boundary responses do not change with further size increases [16].
Gravity was set to zero in this work. The specimen boundaries were defined using smooth ''wall'' elements. The cubic REV may then be loaded in true-triaxial space. Rigid boundaries are smooth so that principal axes of stress and strain are coincident with the cube axes. The principal strains are calculated directly from the wall displacements, while the corresponding principal stresses are obtained from the boundary forces. Any compatible stress-strain path may be followed: strain components may be imposed directly as wall motion, whereas stress components are enforced through servo-control. Ciantia et al. [16] examined the spatial distribution of crushing events during incremental probing and showed that the rigid walls induce no localization.

Stress probing and strain response envelopes
The complete set of simulations carried out in the current study includes some ancillary triaxial tests, but is mostly dedicated to incremental stress probing starting from different initial states on the triaxial plane. For axisymmetric (triaxial) loading, the number of independent stress and strain variables reduces to two, and a convenient graphical representation can be given in a Rendulic plane. Referring to Fig. 2a, in the Rendulic plane of stress increments, the stress probe magnitude is: where r z is the stress in the z direction and r x and r y (the stresses in the x and y direction, respectively) are equal due to the axisymmetric conditions. The stress probe magnitude was established as 1% of the mean principal stress value at the corresponding initial state (from which probing is started). The stress probe direction is defined by the angle a Dr between the ffiffi ffi 2 p Dr x axis and the stress increment vector: At the specimen scale, the response to a set of incremental stress probes is described using an incremental strain response envelope [10,54]. The strain response envelope is plotted in a Rendulic plane for strain increments, (De z : ffiffi ffi 2 p De x ) (Fig. 2b). The strain increment direction is defined by the angle, a De , between the ffiffi ffi 2 p De x axis and the strain increment vector: In Fig. 2b, the strain increment directions corresponding to pure deviatoric strain (De v ¼ 0) and pure volumetric strain (De d ¼ 0) are also indicated.
The resolution at which the strain envelope should be obtained is dictated by the purpose of the study. Ciantia et al. [16] used only 12 probes for each initial condition; this resolution was appropriate for identifying the main behavioural trends but insufficient to accurately describe the material response. In this work, where the fabric as well as the strain responses are analysed, 108 probes were performed at each probing point (36 incremental directions, spaced at 10 deg a Dr , were probed each with 3 parallel probes-see below).

Fabric descriptors
The mechanical coordination number and fabric tensor were systematically computed as descriptors of specimen fabric in this study. The mechanical coordination number is defined for a collection of particles as [57]: where N 1 and N 0 are the number of particles with 1 and 0 contacts, respectively, and N c and N p are the number of contacts and particles in the specimen, respectively. As the definition of Z m cannot be applied to a single particle, when presenting results for single particles we refer to connectivity, defined as the number of contacts per particle. The fabric tensor / ij is given by: where N c is the total number of contacts and n k i denotes the unit contact normal for the k-th contact. Fabric anisotropy is quantified using / d , the second invariant of the deviatoric fabric [4], sometimes referred also as the Von Mises fabric invariant [45]: where / 1 and / 2 are the major and minor principal values of the fabric tensor. Following Thornton and Zhang [58], for a set of incremental stress probes, it is possible to determine an incremental fabric response envelope. However, in this study, the direction of fabric incremental change remains fixed. This is a consequence of the general condition of null fabric trace change D/ 1 þ D/ 2 þ D/ 3 ¼ 0 and that, for axisymmetric loading, the fabric maintains axial symmetry [57] and / 2 ¼ / 3 . In this work, the principal directions 1, 2, 3 coincide with the z, x and y axis directions, respectively, and, for consistency with Eqs. (4)-(6), it follows that ( Fig. 3a) Plotting the results in the equivalent Rendulic plane (D/ z : ffiffi ffi 2 p D/ x ) will not add any further information to this result. For this reason, the micromechanical interpretation of stress probe results is best illustrated by plotting the change in fabric component magnitudes D/ k against the probe stress loading direction a Dr (Fig. 3b).

Strain and fabric response decomposition: parallel probe approach
The so-called parallel probe approach is a numerical technique that applies the same stress probe to initially identical DEM specimens with differently specified contact properties [54,60]. As detailed in Ciantia et al. [16], parallel probes may be used to decompose the incremental strain envelope into a reversible (''elastic'') part, De e , and two irreversible (i.e. ''plastic'') parts: De pu due to sliding between uncrushable particles, and De pc due to particle crushing, so that To achieve this decomposition, three simulations were run for each stress probe. The first simulation used the crushable DEM model described above and the measured strain is De; this is termed the elasto-plastic-crushable (e-pc) probe. In the second simulation, all the mechanisms responsible for plasticity (interparticle sliding and particle crushing) are inhibited to give De e ; this is taken to be the elastic probe. In the final simulation, only crushing is inhibited and the measured strain is De e ? De pu ; this is taken to be the elasto-plastic (e-p) probe.
The potential of the parallel probe approach is exploited further in this work. To begin with, the same type of decomposition may be applied to the incremental change in fabric so that: where k may represent x, y or z. Furthermore, to study the mechanisms associated with crushing, at every crushing event during an e-p-c stress probe, the micromechanical state (fabric, coordination numbers) was stored. As explained below, this allowed systematic comparison with similar quantities recorded at the same situation in parallel e-p probes and, therefore, isolation of crushing-induced fabric effects.

Initial states
The three initial states selected for exploration through stress probes (A, B, C; see Fig. 4a) each belong to a p 0constant triaxial compression stress path performed at p 0 = 52 MPa. These states are characterized by different degrees of stress obliquity (defined as g = q/p 0 , where q is the deviatoric stress), with g = 0 at point A; g = 0.5 at B and g = 1 at C. The triaxial response of the crushable model at this stage is ductile and contractive (Fig. 4b). The same stress path was applied with crushing disabled: the response during shear becomes dilatant and slightly brittle (Fig. 4b).
As a reference, Fig. 4a also includes the grading index I G contours on the triaxial plane, that had been obtained by Ciantia et al. [16]. The grading state index, I G , was introduced by Muir Wood [44] as a measure of grading evolution and is computed as an area ratio of the current grading to a limit grading in the grain size distribution plot (see [16]. All three probing points are located in a region well above isotropic yield (p 0 cr = 40 MPa) (Fig. 4c) in which granular evolution through crushing is already well established. As a reference, Fig. 4c also includes a critical state line (CSL) obtained from several auxiliary triaxial tests, including those used in calibration (Fig. 1). As explained in detail by Ciantia et al. [17], this CSL is that relevant to normally consolidated states such as those considered here. The p 0 -constant stress path accelerates the crushing-induced evolution with respect to that observed along the isotropic path (Fig. 4d). Therefore, more crushing is expected upon incremental probing at point C than at point B, and more at point B than at point A.

Crushing effect on fabric development
Before examining the incremental responses, it is useful to have a general perspective on how fabric evolves in the triaxial plane as crushing proceeds. To do so, we examine the initial values of Z m and / d at the probing points, before incremental loading. To obtain a richer picture, we have extracted relevant data from 12 different initial states attained by following stress paths at fixed stress ratios, g = q/p 0 , of 0, 0.3 0.75 and 1 (other results from these tests were presented in [16]. To highlight the effect of crushing, the same stress states were also obtained in simulations that inhibited crushing. Figure 5a, b indicates on the q:p 0 plane the values of Z m for various stress states of the crushable and uncrushable models, respectively. Comparing these values, the most significant effect observed is that Z m reduces after the onset of crushing. This is clearly visible in Fig. 5c where Z m is plotted against mean stress p 0 normalized by the apparent yield stress on isotropic loading, p 0 cr . For values of p 0 / p 0 cr \ 1, the Z m of the crushable and uncrushable samples is practically identical while for p 0 /p 0 cr [ 1, Z m of the crushable samples is significantly lower than the uncrushable ones. There is also an effect of ductility: at the same stress level the crushable material has suffered larger shear strains and is closer to critical state (Fig. 4b). At the critical state (CS), Ciantia et al. [17] observed a practically unique relation between mean stress level and Z m , that is insensitive to grading. This is again noted here, where the end points (Fig. 5d) are very close for the uncrushable (Z m-= 5.91) and crushable cases (Z m = 5.66). This implies a rapid fall in Z m as the CS is approached for the uncrushable material in a dilatant manner. Indeed, such abrupt falls in Z m for uncrushable particles have been documented before in tests with dilatant responses [31,45]. Figure 6 presents results related to the Von Mises fabric invariant. For both the crushable (Fig. 6a) and uncrushable (6b) cases fabric anisotropy / d increases as the stress obliquity g increases. This is in agreement with previous work [57]. Crushing results in a more gradual development of fabric anisotropy as deviatoric stress increases (Fig. 6c); a similar effect had been noted by Xu et al. [64]. However, as it happened with Z m , the differences reduce as the  As the deviatoric stress level increases, stress probing results in more plastic strain. While the elastic envelope changes little, the e-p-c and e-p strain envelopes increase in size. The uncrushable plastic component is larger than the crushable one at points B and C, but smaller at point A. An isotropic stress path had been followed to point A, and a small amount of crushing takes place along probing directions close to that isotropic loading path (marked IC in Fig. 7a). This memory effect reflecting the previous stress path is also clear in the envelopes at B and C, where the Cumulative distribution function of critical contact force ratios before stress probing, for all particles and for particles that will crush Acta Geotechnica 123 maximum straining is close to, but not coincident with, the previous stress path direction (marked DC in Figs. 7b, c). For both the crushing-induced and crushing-independent mechanisms, the plastic flow direction is independent of the stress probing direction. Plastic flow tends to a critical state type of flow (i.e. De vol = 0) as g increases, with both plastic flow components becoming more deviatoric as g increases; however, it is clear that the plastic-crushable component is always more volumetric than the plasticuncrushable one; that was a general feature also noted by Ciantia et al. [16]. Figure 8 details the orientation and magnitude of plastic flow components for all stress probes at point C (g = 1). In Fig. 8, an g represents the orientation of the plastic flow direction, while an f is the direction of maximum flow magnitude: it is clear that these are not coincident.

Strain response envelopes and plastic flow
Both plastic mechanisms are well described by a cosine shaped curve, although the agreement is better for the crushing-independent plastic mechanism. This is due to the relatively small number of crushing events (46 at maximum) that drive plastic deformation in the crushing-induced mechanism. Fig. 11 CDF curves of the contact normal forces for different contacts subsets: all contacts on crushed particles; FR cr p contacts on crushed particles; all contacts; FR cr p contacts on all particles. For a given particle FR cr p is the maximum of F/F lim , at its contacts (a) (b) (c) Fig. 12 Connectivity CDF curves for different particle subsets: particles that will crush; particles having a FR cr p [ 0.8; all particles. For a given particle FR cr p is the maximum of F/F lim , at its contacts

Micromechanical observations on particles that crush
Particle crushing induces a local dynamic instability in the granular network [16] so that crushing events are typically grouped in time and can be described as ''cascading'' [18,28]. Cascading means that later crushing events are triggered by the earlier ones, through force redistribution.
Using the parallel probe technique, the connected and cascading nature of different crushing events is easily demonstrated. Results from the probe with largest number of crushing events (a Dr = 140 at g = 1) are shown in Fig. 9. For a given particle, we define the critical contact force ratio, FR cr p as the maximum of F/F lim , at its contacts. When crushing is enabled that ratio remains always below 1, but if crushing is disabled it may rise above 1. For instance, Fig. 9a tracks the rise of FR cr p during an e-pc probe on some particles up to their crushing moment, and compares it with the situation of the same particles in a parallel e-p probe, where crushing is disabled. When crushing is disabled, forces on some, but not all, of those particles rise above the crushing limit. The ordering of crushing events is significant in that respect. In Fig. 9b, it is shown that particles that crush earlier in the e-p-c probes are those that experience forces exceeding the crushing limit in the e-p probe, whereas those that crush later in the e-p-c probes remain well below the limit if crushing is disabled.
Particles that will break during probing are those at the top of the FR cr p distribution (Fig. 10). About 75% of particles are very far from breakage, with FR cr p values below 10%. Breakage during probing appears limited to particles that start the probe with FR cr p values above 80%. A high critical contact force ratio is not the same as high contact force. Forces acting on particles that go on to crush are not necessarily those at the top of the contact force distribution. Figure 11 represents the cumulative distribution function (CDF) of contact normal forces for different contact subsets. Forces acting on particles that will crush are more narrowly distributed than forces acting on all particles. The difference increases as shear advances and, as a result of crushing, the number of relatively stronger smaller particles is increased. Several models for crushable soils (e.g. Lobo-Guerrero and Vallejo [41]; Ben-Nun and Einav [5]) have included a connectivity-dependent term in their breakage criteria. In contrast, no link to coordination number was prescribed the breakage criteria implemented here. It was then an open question if particle breakage was in any way related to connectivity. To answer this, particle connectivities for particles that are closest to crushing (i.e. those with FR cr p values above 80%) are compared with that of all particles in Fig. 12. These data indicate that connectivity is not, per se, a good indicator of crushing likelihood. The main difference is that the CDF of particles that are close to crushing appears-again-somewhat more narrowly distributed (i.e. includes less particles with less than 3 or more than 9 contacts) than the general distribution. 2. For a given particle FR cr p is the maximum of F/F lim , at its contacts A complementary illustration of this effect is given in Fig. 13, which illustrates the evolution of connectivity on particles that crush during a representative probe (a Dr-= 140 at g = 0.5, in which the number of crushing events allows a clear plot). Connectivity of particles that crush does not cluster at any particular value and is not much affected by previous crushing events.
Contact fabric data is presented in Fig. 14 for the three initial states before probing. This may be compared with Fig. 15 which is restricted to particle critical contacts, N c,cr (for a given particle, N c,cr is the contact closest to crushing, hence the one with the maximum of F/F lim ratio). As shear progresses force orientation becomes more anisotropic. Critical force contacts (Fig. 15) have a stronger normalized average force magnitude and are more anisotropically oriented that the general contact network. Furthermore, as the shear-induced contact force anisotropy increases, the number of particles whose critical ratio is above a certain threshold also increases (Fig. 16).

Incremental fabric changes due to crushing
The previous section explored the micromechanical conditions that favour crushing. A separate question is: what is the micromechanical effect of crushing, or, using the terminology introduced previously, what are the mechanisms underlying the plastic-crushable irreversible strain? These plastic strains are due to contact slippage and granular network rearrangement. To identify them, a sequential probe technique was developed.
The sequential probe approach is illustrated in Fig. 17. It comprises several steps: 1. From the initial state O an ordinary e-p-c probe (O-B) is performed, to identify the particles that crush during the probe 2. From the initial state O an elasto-plastic probe (crushing disabled) was performed (O-A) and final state (point A) was taken as a reference condition 3. At the reference condition (point A) particle sliding was inhibited and all the particles that would crush in the corresponding e-p-c probe [determined in (i)] where crushed simultaneously, while maintaining the specimen stress constant. Some very small deformations were observed before equilibrium was restored (point A 0 ). These are due to elastic contact force readjustment in the network. 4. Particle sliding (without crushing) was re-enabled. The incremental strain path then moved from A 0 to B 0 . The final strain state obtained (point B 0 ) was close to the ep-c probe final point (point B) showing that this sequential approach gives a strain response very close to the e-p-c response. The advantage, however, is that this post-crushing sliding step allows to identify the contacts where slip is induced solely by particle crushing.
This sequential probe method was systematically applied to explore the nature of crushing-induced particle slippage. In Fig. 18, rose diagrams of the sliding contacts in the uncrushable step (left) and the crushing-induced step (right) are represented for three representative probes. (For instance, the uncrushable step corresponds to path OA in Fig. 17, while the crushing-induced step corresponds to path A 0 -B 0 in Fig. 17). Each bin in the rose diagram is shaded by the normalized force value for that orientation; (forces are normalized by the overall force mean value, f n ). The dashed red line indicates a preferred (or average) direction, a ave , for which sliding is maximum. For the uncrushable case this direction is independent of the stress increment direction a Dr . For contacts that slip due to particle crushing that preferential direction is not so obvious and sliding contacts are more isotropically distributed.
The incremental fabric changes due to crushing-induced sliding and to uncrushable sliding are now compared. Results for the strong contact subset network (i.e. network of contacts carrying above-average forces) are presented for both steps of sequential probes in Fig. 19. Restriction to the strong subset network was necessary as incremental fabric changes for the weak contact network were imperceptible (see Sufian et al. [52] for a more general discussion of fractional networks). Figure 19 hence compares the change in vertical strain with the change in the vertical probing direction. For instance, when g = 0 at point A, D/ z = 0 for probing directions (a Dr ) close to IC and IE, while D/ z is maximum/minimum for a probing direction close to the DC/DE, respectively. Probing from the isotropic state resulted in very small plastic deformations (Fig. 7a). This is also reflected in Fig. 19a, where the changes in vertical fabric are almost unaffected by activation of the plastic mechanisms. The situation at point B (g = 0.5) is different (Fig. 19b): plastic strains are already significant and are accompanied by perceptible fabric changes. The dominant component is non-crushed induced plastic sliding, which also appears to fully control fabric changes. Although crushing-induced plastic strains are not negligible (see also Fig. 7b) they seem to leave fabric unchanged. When g = 1, at point C (Fig. 19c), the trend towards increased crushing-induced plastic strains is more clear. The effect on fabric of those crushing-induced strains is to reduce the changes that pure frictional sliding will induce, i.e. to reduce the development of anisotropic fabric.
This effect is clearly visible in Fig. 20 which considers only the incremental fabric changes due to plastic-uncrushable strain (D/ pu z ) and plastic-crushable strain (D/ pu z ) for probes starting at g = 1. Each incremental fabric change response is oriented almost identically to the corresponding incremental strains: compare the maximum strain direction in Fig. 8 and that of maximum incremental fabric change in Fig. 20. This directional coincidence is, perhaps, less surprising than the different sign of the fabric change components. It is clear that for the purely frictional, uncrushable mechanism, incremental fabric changes are positive, whereas crushing-induced plastic structural mechanism the change in vertical fabric is negative. This means that, since D/ z ¼ À 2D/ x , crushing induces the creation of new strong contacts in the direction perpendicular to the direction of loading and the disappearance of strong contacts parallel to the direction of loading.

Conclusions
This work has examined several aspects of the microscale response of a crushable soil discrete element model. The most significant findings are 1. When shearing is accompanied by particle crushing, the evolution of microstructural measures like the mechanical coordination number or Von Mises fabric towards their critical state value is smoothed. A ductile specimen-scale response is accompanied by a slow micromechanical evolution. 2. Plastic deformation due to breakage takes place as breakage-induced contact sliding. Sliding due to this mechanism is more isotropically distributed than noncrushable sliding and tends to reduce fabric anisotropy instead of increasing it. The more volumetric nature of crushing-induced plastic flow is thus explained. 3. In this model particles that break are those in which one contact force overcomes a breakage limit that depends both on particle and contact properties. There is not a clear connection between particle connectivity and particle breakage. On the other hand, increased b Fig. 19 e, e-p and e-p-c vertical strain (left) and vertical fabric change (right) as a function of probe loading direction a Dr at the three probing points. The numbers next to the e-p-c curves represent the number of crushing events fabric anisotropy does take more particles closer to breakage. Shear-enhanced breakage appears then as a consequence of shear-enhanced fabric anisotropy.
Although the third point is clearly dependent on the specification of the breakage model, the second and first points noted above are not. Indeed, breakage smoothing of fabric evolution had been noted in very different discrete breakage models such as those based on aggregates. As the capabilities of experimental micromechanics to address breakable soils improve some or all of the micromechanics uncovered here may be directly verified. In the meantime, performance on mesoscopic simulation would remain the best guide to judge the fitness of discrete models. The role of micromechanics studies, such as this one, is to elucidate the causes underlying a good or a poor performance of the model.