Micro mechanics of drained and undrained shearing of compacted and overconsolidated crushable sand

It is well known that particle crushing has a major influence on the physical behaviour of granular soil. For a sand being sheared, the crushability of the particles affects both the strength and dilatancy of the soil (Hardin, 1985). There have been several efforts to use the discrete-element method (DEM) to model and investigate the general shearing behaviour of soil using crushable particles (e.g. Abe & Mair, 2005; Lobo-Guerrero et al., 2006; Wood & Maeda, 2008; Marketos & Bolton, 2009; Elghezal et al., 2013; Wang & Yan, 2013; Ciantia et al., 2016). In particular, a number of recent attempts have been made to investigate the three-dimensional triaxial behaviour of crushable soil (e.g. Bolton et al., 2008; de Bono & McDowell, 2014; Hanley et al., 2015; McDowell et al., 2015). However, in general the majority of published studies of soil behaviour using DEM neglect particle crushing (Thornton, 2000; Sitharam et al., 2002; Cui & O'Sullivan, 2006; Minh & Cheng, 2013).

One method to take into account particle crushing is the use of agglomerates (e.g. McDowell & Harireche, 2002), which involves representing single grains by bonded groups of spheres. Discrete-element modelling of triaxial tests using such agglomerates (Cheng et al., 2004; Bolton et al., 2008) provides useful quantitative insight on aspects such as yielding and plastic deformation, as well the general effects of particle crushing during shearing. However, agglomerates are problematic for a number of reasons, principally due to the finite limit on fragment size and difficulty with measuring/interpreting voids ratio. This can be overcome by using the alternative approach of replacing breaking particles with smaller particles of the same shape, such as by de Bono & McDowell (2014) and Hanley et al. (2015), who both provided overviews of the effects of crushing, and established critical state lines (CSLs) from triaxial tests. Ciantia et al. (2016) also performed shear tests using the particle replacement method; however, their study focused on the plastic flow and response to strain probes, rather than critical states.

The authors have published several recent studies using a simple crushing model, and have previously attempted to establish a CSL (de Bono & McDowell, 2014; McDowell et al., 2015). This paper follows on from that work by aiming to establish a full CSL over a wide range of stresses for a simulated silica sand and provide an in-depth fundamental analysis of the micro-scale behaviour of a crushable soil during and after being sheared to a critical state, from initial states looser and denser than critical, involving a range of different stress paths, including constant-volume conditions.

The simulations presented here are performed using the software PFC v5 (Itasca, 2015). The results are all obtained from a single initial sample of spherical particles, shown in Fig. 1. The initial cylindrical sample is subjected to isotropic normal compression to high stresses. Triaxial shear tests are performed from various states of isotropic compression, as well as from a range of overconsolidated (OC) states.

This work extends the use of the model proposed by McDowell & de Bono (2013) for normal compression to the study of drained and undrained stress paths on compacted and OC samples and the establishment of a state boundary surface. To permit the extensive number of simulations required in an acceptable time frame, the initial sample contains a relatively small number of particles – 723. The choice of such a small quantity of initial particles is due to the extensive crushing that occurs in nearly all of the simulations, resulting in far greater quantities over a range of sizes. The simulations presented here have no comminution limit – that is, there is no lower limit to particle size. This approach allows the unadulterated evolution of a fractal particle size distribution (PSD), and for particles to continue getting smaller. This means, however, that the numerical time step used in the simulations – which is dominated by the square root of the mass of the smallest particle – becomes very small. Hence, to enable the number of simulations required here while not imposing artificial conditions, the initial number of particles is small, which helps to minimise the calculation time (i.e. the time taken to compute the motion, and so on, for all particles) once extensive crushing has occurred.

Other researchers often impose an arbitrary comminution limit with no physical basis in order to manage computational time. For example Hanley et al. (2015) started with approximately 22 000 initial particles, with an apparent size range of 0·13–0·62 mm. They imposed a lower limit to particle size of 0·12 mm, fractionally lower than the initial smallest size. Using a high-performance computing facility, their simulations at the highest pressures with the most crushing resulted in around 40 000 particles.

In comparison, the most extreme case presented here consists of around 17 000 particles after shearing, and has d_max/d_min = 40. These particles are generated from the original 723. Repeating this simulation using a sample of, say, 20 000 initial particles would lead to many hundreds of thousands, if not more, of very fine particles, which is simply not computationally possible.

The initial sample has a ratio of particle diameter to sample diameter (l/d_max) of 15/2. Although this is unattractive, from an experimental point of view, in particular with regard to testing of rockfill and railway ballast, it is commonly accepted that a ratio of l/d_max of 6/1 is enough to eliminate any sample size effects on the measured results (e.g. Marachi et al., 1972; Indraratna et al., 1993; Lackenby et al., 2007). That said, it should be noted that, in nearly all of the simulations presented here, crushing increases the number of (smaller) particles greatly. With the exception of the simulations at the lowest pressures which do not undergo crushing (and will be highlighted), at critical states the samples have total particles numbering 2250–16 500. Using a general descriptor of particle size such as d₁₀ or d_ave (e.g. as suggested by Hardin (1989) and by Marketos & Bolton (2010)), the simulations give values in the ranges 0·47 < d₁₀< 0·95 and 0·52 < d_ave< 1·25 at critical states. As such, the sample used here may be considered as representative, and therefore there is no need for periodic boundaries – the material here is compared with established, real laboratory data, measured at the boundaries, and for this purpose a flexible membrane is used. In addition, future work will focus on the effect of particle shape, which should lead to shear bands for heavily OC sands in the present authors’ simulations; these would be prohibited by periodic boundaries and, to be consistent, are not used here.

Furthermore, the normal compression line (NCL) presented in this work is compared to one obtained previously using a sample with a different aspect ratio and a higher percentage of boundary particles in a rigid oedometer (de Bono & McDowell, 2016a), and there are no differences to the shape or slope of the NCL. Also, the ability to create and reuse any specific numerical sample in the DEM mitigates any concerns to do with representative samples, and means that a DEM sample containing as many particles as a standard laboratory sample is not required.

For computational efficiency all particles are modelled using spheres, and gravity is neglected. Key parameters are given in Table 1. The sample is enclosed vertically by two rigid horizontal platens, and laterally by a rigid but deformable ‘cylindrical’ boundary. This lateral ‘membrane’ is constructed from triangular facets, whereby the vertices can move independently from one another in order to ensure a uniform pressure is applied to the specimen. The number of facets affects the computational time required to ensure uniform application of confining pressure during the simulation; in this work a total of 4320 facets was deemed a suitable compromise between membrane flexibility, confining pressure homogeneity and computational time.

Axial strain (ε_a) is applied by gradually accelerating the upper platen to a velocity of 0·01 m/s. The membrane vertices are also given a vertical component of velocity proportional to the platen velocity (i.e. using a constant strain rate). The axial stress σ_a is measured as the sum of all normal contact forces acting on the loading platen divided by its area (averaged for both top and bottom platens). The radial stress σ_r is equal to the applied confining pressure.

The sample volume (and associated strain) is calculated by dividing the space enclosed by the walls into a series of pyramids, where each membrane facet (and both platens) is the base of a pyramid which has its apex in the centre of the specimen. The volume of these constituent pyramids can be easily calculated and summed. The axial strain is calculated directly from the position of the upper platen, and together with the volumetric strain is used to estimate the current radial strain ε_r, assuming the sample approximately remains cylindrical in shape (the common assumption). The triaxial shear strain is calculated as

Simulations are typically run until either an axial strain of 30% or a critical state is achieved, with the exceptions of some of the computationally demanding simulations (e.g. those at the highest stresses with substantial crushing).

The authors’ previous work is based on the same crushing model, which replaces broken particles with new smaller fragments while obeying conservation of mass. Key aspects of this particle breakage mechanism will be summarised as follows. Initially, all particles are given randomly distributed strengths, in terms of octahedral shear stress, q, measured as

where σ_1,2,3 are the average principal stresses in the particle, caused by and calculated from contacts with neighbouring bodies. The particle strengths for a given size are distributed according to a Weibull distribution, defined by a characteristic value q₀ and a modulus m, which determine the location (approximate average strength) and variability of the distribution. The characteristic value q₀ is a function of particle size according to

where b represents the size effect on strength, and in this case, using Weibull statistics

The characteristic strength q₀ is a value such that 37% of particles are stronger, and is useful as a gauge of the average strength for a particular particle size (it is similar to the mean value of the distribution). The particle strengths and Weibull modulus (Table 1) are the same as used previously by the authors and are obtained from experimental particle crushing tests (McDowell, 2002).

The breakage model is identical to that used in the authors’ previous work, to which readers are directed for full details (McDowell & de Bono, 2013), and which includes discussion on the number of fragments and hardening law; readers are also referred to more recent work, which provides an in-depth study on the choice of breakage criterion (de Bono & McDowell, 2016b).

For application of an isotropic pressure, the simulation sequence starts with applying a small stress increment using the wall servo-controls. Once an increment is completed (and the walls are stationary), all particle stresses are checked, and any particle in which the octahedral shear stress exceeds the strength is replaced by two smaller fragments. These fragments are placed within the boundary of the ‘parent’ particle, hence the new fragments overlap. They are aligned in the direction of the minor principal stress (e.g. McDowell & de Bono, 2013). Next, a number of time steps are completed which allow overlapping fragments to move apart, until the sample reaches stability and all artificial overlaps have been dissipated. If particle breakage results in a drop in the applied isotropic stress then the stress increment is reapplied. Once an isotropic stress is sustainable without particles breaking, the process is continued (i.e. the next stress increment is applied) until the desired macroscopic stress state is achieved.

For the triaxial shear tests, a similar procedure is followed, the notable difference being that these simulations are strain-controlled. After a particular isotropic stress has been applied to the sample, the conventional triaxial shear tests begin by applying an increment of axial strain (0·1%) to the sample. Simultaneously, the membrane servo-control is activated to ensure that the confining pressure σ_c remains constant, alleviating any increase in radial stress that would be caused by the imposed strain. After the axial strain increment and the correct uniform radial confining pressure are applied, all particles are checked and allowed to break, as outlined above. After any particles have broken, the radial confining pressure is reapplied if necessary. Once breakage (if any) is complete and the correct confining pressure is applied, the process continues to the next strain increment.

In addition to conventional triaxial tests, constant-p′ and constant-volume triaxial tests are also performed. For the constant-p′ tests, the simulation procedure is similar to that outlined for the standard triaxial tests, but the membrane servo-control achieves a specific mean stress (p′) rather than confining pressure. For the constant-volume simulations, likewise the radial stress is adjusted to achieve the correct sample volume.

The results from the isotropic normal compression simulation are presented in Fig. 2. Included in this figure is the NCL obtained from one-dimensional normal compression of the same material (de Bono & McDowell, 2016a). The two NCLs exhibit the same slope, and are parallel, with the one-dimensional NCL to the left of the isotropic NCL, as is observed experimentally; although the separation is evidently small in the simulations due to the lack of particle shape, giving limited ability for particle interlock and leading to a high K₀ value. The slope of the isotropic normal compression simulation is approximately the same as the one-dimensional test, and gives a slope of −0·5 on log–log axes, consistent with previous work and predicted by the authors’ compression law (McDowell & de Bono, 2013)

where e and p′ are the current voids ratio and stress, respectively, e_y and $p′y$ are the voids ratio and stress at yield, and b is the size effect on particle strengths from equation (3) – that is, b = 3/m. The term 1/2b describes the slope of the compression line, and using a value of m = 3·3 for silica sand (McDowell, 2002), this predicts a slope ≈ − 0·5, indicated by the dashed line in Fig. 2(a). Fig. 2(b) shows the points from where triaxial tests are performed. Triaxial tests performed from points that are before yield (<15 MPa) (and which may lie on any non-unique ‘virgin’ loading line) are referred to as ‘compacted’, following the terminology of Coop (e.g. Coop, 1990; Jovičić & Coop, 1997). This is to distinguish from states that are beyond yield (≥15 MPa), on the high-stress, unique NCL, which are referred to as normally consolidated (NC).

Conventional triaxial tests were performed from isotropic stress states of 1, 3, 5, 8, 10 and 15 MPa; the results are given in Fig. 3. Only the test at 15 MPa is from an NC state. No conventional triaxial tests at higher pressures were conducted because of computational limitations. At large strains at such high pressures, particle crushing results in a highly developed PSD with many small particles, meaning such simulations were very time consuming to continue. Images of the sample from the conventional σ_c = 15 MPa test are shown in Fig. 4, at ε_a = 0 and 35%. At the end of the test, the sample shows extensive crushing, with approximately 14 000 particles (d_min= 0·05 mm) and barrelling failure (due to the ability of the membrane to accommodate deformation).

The scatter evident in these plots is due to the particle breakage mechanism, which causes an immediate loss of contacts in the major stress direction. This could potentially be overcome by using vastly greater quantities of particles, but this is not currently possible – even the samples which experience the most breakage, and contain many thousands of particles, still demonstrate such fluctuations (e.g. the σ_c = 15 MPa test).

In Fig. 3, all simulations reached states of approximately constant stress, with mild peak states evident at the lowest pressures and purely ductile behaviour at the highest. Increasing the confining pressure suppressed dilatancy, and the tests at the highest pressure required large strains to approach a state of approximately constant stress.

Constant-p′ tests were performed from 1, 3, 5, 8, 10, 15 and 20 MPa, and similar trends were observed. The results in general compare well to typical experimental data (e.g. Yamamuro & Lade, 1996).

Constant-volume tests were performed from initial states of 8, 10, 15, 20 and 25 MPa; the results are given in Fig. 5. For plotting the data in Fig. 5(b), the stress responses have been smoothed in order to clearly see the intrinsic stress paths. Plotting the raw data in this form obscures the actual stress path due to fluctuations about the critical state. The smoothing technique used is similar to a ‘moving average’ technique, and is used solely to allow clear q–p′ plots; it does not in any way affect any of the analysis presented here.

The simulations are consistent with experimental data, and demonstrate typical undrained characteristics such as phase transformation points (PTPs) (Ishihara, 1993). A PTP is a feature of medium-to-dense undrained experimental tests (e.g. Been et al., 1991; Ishihara, 1993) and represents the stage at which the material behaviour changes from contractive to dilative – that is, when the rate of change of p′ reverses.

In the compacted test at the lowest pressure, the deviatoric stress increases monotonically until reaching a stable value at around ε_a= 6%. The mean stress meanwhile initially decreases until reaching the PTP at about ε_a= 2%, after which it increases, together with q, at an approximately constant stress ratio; eventually it reaches a critical state. It was not possible to conduct ‘undrained’ simulations at lower confining pressures, due to the initial compressive response at the start of shearing – there was not enough scope for the confining pressure to reduce sufficiently to maintain the original volume. The servo-control in such tests causes the stress to drop to 0 in pursuit of the correct volume. In other words, such simulations display rapid dry liquefaction with strain softening, as observed in experimental data (e.g. Been et al., 1991).

Increasing the confining pressure (i.e. the initial isotropic stress) suppresses the dilatant tendency of the material; the (NC) test performed at the highest pressure (25 MPa) displays a large decrease in mean stress as the state approaches a critical state. This simulation displays no clear PTP or dilative behaviour, consistent with experiment data for contractile tests, although typically one would also expect an instability point – a local maximum in q in the q–p′ plot (Ishihara, 1993) – behaviour which is not clearly evident. In general, the transition from dilatant to contractile behaviour appears to occur at a relatively higher stress level in the simulations compared to most available experimental data, and as such, it is expected that simulations at higher stresses would be required to observe such instability points. It is worth noting that the rapid drop in mean stress p′ evident for the highest pressure simulation (25 MPa) is due to a statistically large amount of breakage that occurs immediately upon commencement of shearing. Such contractive behaviour is only possible when using crushable particles, using unbreakable spheres would consistently result in dilatant responses (e.g. Sitharam & Vinod, 2008).

The isotropically compressed sample was unloaded from 20 MPa, and then conventional, constant-p′ and constant-volume OC tests were performed from stresses of 5, 8, 10 and 15 MPa along this unloading line. In addition, two further high-stress tests were performed after unloading from 30 MPa, at stresses of 15 and 20 MPa. The unloading lines and states from which these OC tests are performed are shown in Fig. 2(b). The results from these tests were more dilatant overall; for example, the 5 MPa conventional simulation displayed more volumetric dilation than the compacted test at the same pressure, and the 15 MPa test displayed less contraction than the equivalent NC test – although it still displayed a large volume reduction.

The compacted, NC and OC tests all reached ultimate stress ratios of around 0·7, with the exception of some of the compacted samples at the lowest stresses, in which no crushing occurred, which gave lower ultimate η-values of around 0·6. Fig. 6 shows the ultimate states – referred to as critical states from here onwards – in q–p′ space, for all simulations. Also shown are the peak states from the conventional and constant-p′ tests, and the PTPs from the constant-volume tests. A CSL with a constant M = 0·68 appears to fit all critical states well. The PTPs indicate the transition from contractive to dilative behaviour, and are often assumed to lie on the CSL, by the same reasoning that points of minimum volume in drained tests correspond to the critical state stress ratio (e.g. Coop, 1990). Taking the PTPs as the minimum values of p′, these points also appear to agree well with the estimated CSL. The peak states display a curved failure envelope, consistent with experimental (e.g. Bolton, 1986) and numerical results with crushable particles (Cheng et al., 2004). The peak envelope appears larger for the OC tests in Fig. 6, indicative of increased dilation. For both sets of data, this envelope appears to converge to the CSL at high stresses, for the tests that are near or on the NCL.

Figure 7 shows the critical states in log e–log p′ space. Fig. 7(a) shows the critical states from the compacted and NC tests, and the first observation is that the CSL appears ‘bilinear’ or curved, and has the same shape as the NCL. The transition point in this case occurs around p′ = 10 MPa, and, similarly to how the yield stress on an NCL indicates the onset of particle crushing, this stress is the point at which crushing becomes the dominant factor in determining the volume of the soil. When sheared well below this stress, there is no crushing and particles are able to sustain inter-particle contact forces as they shear over and past one another as the sample deforms. Beyond this stress, inter-particle contact forces are high enough to induce particle breakage, leading to contraction and state paths in e–p′ space that travel downwards to lower voids ratios (Wood, 1990; Been et al., 1991). The second key observation is that the CSL is parallel to the NCL at high pressures, as is observed experimentally (e.g. Schofield & Wroth, 1968; Coop, 1990), and therefore has a slope of approximately −0·5, indicated in the figure.

Figure 7(b) shows the critical states for all tests, including those from OC states. The critical states from the OC tests also form a curved line with a clear change in slope. The OC triaxial tests at high pressures appear to reach the same CSL as the high-pressure NC tests. The OC tests at low pressures, however, appear to define a new CSL at a lower voids ratio than that for the compacted samples at low stresses. In other words, the compacted and OC samples produce two different CSLs that merge at high stresses. This is the same as observed both experimentally and numerically (Wood & Maeda, 2008; Bandini & Coop, 2011), and is due to the change in grading caused by crushing. Comparing compacted and OC samples at the same stress – the distance between the initial and critical states are greater for the OC samples than for the compacted samples, reaffirming the increased dilation observed in the stress–strain plots.

The elastic unloading and loading behaviour in the simulations can be described fairly well by a constant elastic bulk modulus. Assuming a constant (measured) elastic modulus of 900 MPa on unloading (shown in Fig. 8), together with a Poisson ratio of 0·24 obtained from conventional triaxial tests, assuming isotropic elasticity, gives a shear modulus

of 566 MPa. These parameters were used to calculate elastic strains at each point during shearing, which were subtracted from the total strains to provide the plastic strains. The plastic strain increment vector $δεvp/δεqp$ was then obtained throughout all tests as the gradient from plots of plastic volumetric strain against plastic shear strain.

Figure 9 shows standard stress–dilatancy plots (e.g. Coop, 1990; Wood, 1990) in which the stress ratio is plotted against the ratio of plastic strain increments (⁠$δεvp/δεqp$⁠). Each plot shows results for compacted (dilatant) and NC (contractile) tests for each type of test, in which the aim is to show the effect of increasing stress level. Zero volume change occurs at a stress ratio of about 0·7 for all simulations, which is also the final stress ratio in all simulations. In Fig. 9(a), increasing confining pressure causes the maximum dilatancy (negative) to reduce and the initial plastic volumetric contraction at low stress ratios to increase, as the initial/confining stress increases from 8 MPa (pre-yield) to 15 MPa (post-yield on the NCL). Similar behaviour is evident for tests under increasing confining pressures under different stress paths (Figs 9(b) and 9(c)), even for the constant-volume tests (Fig. 9(c)) in which the elastic and plastic strains are equal; this lends credibility to the method adopted here of calculating plastic strains. At high stresses the stress–dilatancy behaviour becomes approximately independent of confining stress, suggesting that these NC samples are following the state boundary surface; this will be explored further shortly.

Figure 10(a) shows the results for the conventional shear tests at 15 MPa but with different stress histories (p₀ denoting the preconsolidation pressure). The heavily OC soil exhibits clear peak strength and dilation. The peak strengths and strain increment ratios for all simulations are plotted in Fig. 10(b). In this figure the most commonly used stress–dilatancy equations are also plotted: Cam Clay, modified Cam Clay and Rowe (1962). The most notable outliers (indicated) are those simulations at the lowest stress levels which exhibited no crushing. All three relationships overestimate peak strength for most dilatant tests, but this is a well-known feature of these models. The ability of these continuum models to capture the general behaviour of the crushable DEM soil is pleasing, and in particular, comparing Figs 10(a) and 10(b), it can be seen that a conventional triaxial test on the NC soil at 15 MPa would follow the general path given by these continuum models. However, a more general equation of the form

as used by Coop (1990), may describe the data more closely, and this is shown in Fig. 10(b).

In order to plot the state boundary surface, it is necessary to normalise the deviatoric and mean effective stresses by the isotropic preconsolidation pressure, p₀, throughout all tests. This is calculated according to Fig. 8 by plotting the unload–reload line for the current state in question (point A) and extrapolating back onto the isotropic NCL. For the compacted soils (those on the virgin line before yield) p₀ is calculated using the same method, assuming a linear isotropic NCL in log e–log p′ space; that is to say the linear part of the NCL is extrapolated to low stress levels.

The normalised stress paths are plotted in Fig. 11(a) for all simulations and the results are quite remarkable. The DEM simulations using this very simple crushable model have generated a well-defined state boundary surface, which show (and limit) all possible combinations of stress that are possible to be reached with this soil. A typical result for crushable soil by Coop (1990) is plotted in Fig. 11(b) and the ability of the DEM model to reproduce similar behaviour is pleasing. The shape of the state boundary surface suggests a simple two-mode behaviour, after Wood (1990), with a simple shearing Mohr–Coulomb type mode of behaviour at high stress ratios and a volumetric cap at low stress ratios. The strain increment ratios were then used to plot the flow vectors on the same axes as the associated stresses in Fig. 12, alongside comparable experimental data for sand from Coop (1990) (note that this was for a carbonate sand, so the preconsolidation pressure is much lower). The DEM crushable soil clearly has a non-associative flow rule, in agreement with available data for sand, and at least for these data, the critical state does not quite appear to be at the normalised peak deviatoric stress, again in agreement with available data (e.g. Chandler, 1985; Wood, 1990).

The stress paths for two constant-volume tests sheared from 15 MPa with different stress histories are shown in Fig. 13 alongside comparable experimental data from Coop (1990) for carbonate sand. In both cases, the OC sample's stress path rises initially much more steeply, whereas the compacted samples display a rapid decrease in mean effective stress due to the increased particle crushing that occurs.

In order to add further credibility to the approach, the normalised stress paths are plotted in Fig. 14 for an NC sample with an initial mean effective stress of 15 MPa, for the three alternative stress paths. The normalised plots have been deduced by calculating p₀ throughout all tests according to Fig. 11(a). The end states for all three tests lie on the linear part of the CSL and hence the normalised stress paths all reach the same point in Fig. 14.

In contrast, for the OC tests in Fig. 15, the soils compressed to 20 MPa then unloaded to 8 MPa were sheared to critical states which lie on a non-unique CSL and thus the normalised stress paths do not coincide.

Finally, the normalised stress paths for three conventional tests at 15 MPa are shown in Fig. 16, for initial states which are NC and OC with two alternative preconsolidation stresses of 20 MPa and 30 MPa. The NC test follows the state boundary surface from the beginning; the heavily OC soil exhibits mainly elastic behaviour before reaching/following the state boundary surface to a critical state; and the lightly OC soil displays intermediate behaviour. The ability of the crushing model to capture the essential features of soil plasticity shows promise.

It is known that the decrease in voids ratio as a sample moves down the NCL is a direct result of particle crushing, with a wider distribution of sizes and smaller particles required for the soil to be able to occupy a denser state. As such, it can be expected that similar principles apply on the CSL, and that a point on the CSL should have a unique PSD no matter what path is taken to reach it. This appears to be confirmed by Fig. 17(a), which shows the PSDs from three tests with different stress paths that reach similar critical state voids ratios. The tests shown – conventional at 15 MPa, constant-p′ at 20 MPa and constant volume at 25 MPa – have final voids ratios of 0·508, 0·503 and 0·526, respectively; their paths are indicated in Fig. 7(b). The constant-volume test appears to have a less-developed PSD, which correlates with the higher voids ratio and supports the idea that the critical state voids ratio is a function of the PSD and independent of stress path. Also shown is the PSD from a state on the NCL at e = 0·507, which shows less crushing than the sheared samples at approximately the same e.

Figure 17(b) displays the final PSDs for all triaxial shear tests on a conventional plot and Fig. 17(c) shows in terms of number on log–log axes. The former shows that, in general, crushing in all simulations has a similar result on the PSD, with increasing stresses producing the most developed grading. Fig. 17(c) shows the same distributions, but in terms of quantity, on logarithmic axes to assess the fractal character. It appears that, just like normal compression (e.g. McDowell & Daniell, 2001), crushing during shearing leads to a fractal distribution with a fractal dimension of around 2·5, in this plot indicated by the slope.

Strictly speaking a fractal distribution would be infinite (i.e. span all particle sizes), but in reality a fractal PSD for a soil must be limited (McDowell & Bolton, 1998), and is bounded by a largest and smallest particle size. The largest particle size in all cases here is 2 mm, whereas the smallest fractal size depends on the extent of crushing. Hence, in Fig. 17, the smallest sizes for each of the simulations are still emerging, which is a gradual process, and do not adhere to the ideal fractal shown by the trend line. Those simulations at the highest pressures undergo the greatest crushing and have the broadest fractal PSD.

Insight can be gained by monitoring the average coordination number, which in this work is measured as the mechanical average coordination number (Thornton, 2000). This is the average number of contacts per particle, discounting mechanically redundant particles with zero contacts (due to the lack of gravity).

It was found that all simulations approach an average coordination number of approximately 5·1–5·2 at a critical state; examples are given in Fig. 18 for compacted dilatant samples (Fig. 18(a)) and NC contractile samples (Fig. 18(b)). In similar studies, unique values of critical state coordination numbers have also been observed (e.g. Thornton, 2000; Sitharam & Vinod, 2008; Salot et al., 2009), although in some studies (e.g. Wood & Maeda, 2008; Maeda et al., 2010) the critical state coordination number was shown to increase with stress, in contrast to the observation here, which suggests a unique value independent of stress level. However, none of the aforementioned studies considered particle crushing.

Figure 18(a) shows the results of three different tests from a single compacted state, and all three stress paths reach a similar ultimate coordination number, in the region 5·1–5·2. What is interesting here is that the constant-volume example shows a rapid initial decrease from 5·3 to about 5·1, which coincides with the PTP at the critical stress ratio (minimum p′). At this point the mean stress p′ begins to increase, and the coordination number then peaks slightly before approaching a value of around 5·2.

Figure 18(b) shows behaviour for the three types of simulations from an NC state of 15 MPa. On the NCL after yield, the sample typically has a high average coordination number, which reduces upon shearing, even when a sample contracts. This plot shows that, from the same NC state, shearing by way of different stress paths has an almost identical effect on the coordination number, reducing from about 5·6 to around 5·1, despite the final voids ratios being different in these particular tests. This suggests that there is a unique value of average coordination number at a critical state, independent of stress level.

The average coordination number is plotted as a function of p′ in Fig. 19(a) for both normal compression and critical states. For isotropic normal compression, there is a sharp rise up until yielding occurs, thereafter the average coordination number drops and remains approximately constant at a value of 5·6. During the initial elastic behaviour, the coordination number appears to correlate with voids ratio and/or stress, similarly to previous observations (e.g. Sitharam et al., 2002; Wood & Maeda, 2008). A similar trend can be observed for the critical states: some variation at low stresses, before particle crushing becomes prominent, thereafter an approximately constant value of 5·1–5·2 is visible, which is apparently independent of stress level (and also therefore the number of particles/PSD) and would appear to be related to the macroscopic stress ratio.

During normal compression, the strength of the smallest (strongest) particles is proportional to the applied stress (de Bono & McDowell, 2016a). Recalling the hardening law in equation (3) therefore provides a link between applied stress and the size of the smallest particles

For a fractal distribution with D = 2·5, the number of particles (of size L) equal to or greater than a size d_i can be expressed as

and hence using the smallest particle size, d_s, gives an expression for the total number of particles, N_t

Substituting equation (8) into equation (10) therefore gives an expression for the total number of particles as a function of applied stress (de Bono & McDowell, 2016a)

To test this, the total number of particles is plotted as a function of mean stress for both isotropic compression and critical states in Fig. 19(b), and it can be seen that this relationship appears correct for both critical and NC states (at high stresses), with noticeably more particles for a critical state compared to NC states at the same mean stress, due to the difference in stress ratio causing much greater particle stresses and accelerated crushing.

A cylindrical sample of spheres was isotropically compressed to high stresses using a new model, which features a flexible membrane allowing triaxial tests to then be performed. The isotropic NCL established using this model was parallel with the NCL obtained from previous oedometer simulations using a rigid-walled sample.

Conventional, constant-p′ and constant-volume triaxial tests were performed from an array of compacted, NC and OC stress states in order to establish a CSL. The constant-volume simulations displayed behaviour typical of experimental undrained triaxial tests, and highly contractile results were obtained at high stresses due to the inclusion of particle crushing. The critical states for all simulations (and PTPs for constant-volume tests) gave a unique critical state constant M.

In log e–log p′ space, all triaxial tests approached a unique CSL at high stresses, parallel with the NCL. At low stresses, non-uniqueness was observed, with the OC samples revealing a new, lower CSL, consistent with experimental findings.

In addition, the simple model was able to capture the essential features of the stress–dilatancy behaviour of crushable soil, and was found to give results which are consistent with continuum theories. In particular, a state boundary surface has been exposed, which resembles that for real crushable soil. This lends credibility to the model and the normalised stress paths were therefore plotted for NC and compacted soil, and for OC soil.

Analysis of the PSDs after shearing showed that, similarly to normal compression, shearing leads to fractal PSDs with a fractal dimension of approximately 2·5.

This paper has shown that the simple McDowell & de Bono (2013) crushing model can exhibit all the salient features of soil in states looser and denser than critical, whether compacted or OC and whether drained or undrained. The ability of this model to capture many features of plasticity holds much promise in application to boundary value problems such as pile driving.

This work was supported by the Engineering and Physical Sciences Research Council (grant number EP/L019779/1).

b

size effect on particle strength

d

particle diameter

d_s

smallest particle size

e

voids ratio

G

shear modulus

K

bulk modulus

M

critical state constant

N

number of particles

N_t

total number of particles

p′

mean effective stress

p₀

preconsolidation pressure

q

octahedral/deviatoric shear stress

q₀

characteristic (average) strength for particles of given size

δε^p_p

plastic volumetric strain increment

δε^p_q

plastic shear strain increment

ε_a

axial strain

ε_q

shear strain

ε_r

radial strain

ε_v

volumetric strain

η

stress ratio (q/p′)

ν

Poisson ratio

σ_a

axial stress

σ_r

radial stress