Micro mechanics of the critical state line at high stresses

A critical state line is presented for a crushable numerical soil, which is parallel to the isotropic normal compression line. A previous theory for the normal compression line, which correctly predicts the slope as a function of the size-effect on particle strength is extended to justify the slope of the critical state line. The micro mechanics behind critical states are examined, leading to a theory for a relationship between the volume of smallest particles and mean effective stress. A unique relationship exists for crushed states, leading to a two-dimensional interpretation of the state boundary surface for soils looser than critical.


Introduction
Following Pestana and Whittle [1], McDowell [2] proposed analytically that the high-stress normal compression line for a sand should be linear on log e-log σ axes. This was based on the kinematics of particle crushing and the assumption of a fractal particle size distribution. This was investigated numerically using the discrete element method [3], and it was found that the normal compression line (NCL) for a sand could be described by: where e is the current voids ratio, σ′ is the current stress, b describes the hardening rate of the particles, and e y and σ′ y are constants defining the yield point. The parameter b defines the rate at which average particle strength increases with decreasing size: and as shown by Eq. (1), b determines the slope (=1/2b) of the NCL on logarithmic axes. This compression law in Eq. (1) was tested against experimental data for which normal compression lines and particle strength data (i.e. values of b) obtained from particle crushing tests were readily available in [4], and it was found to correctly predicted the slope of the NCL for each case. It is a well-known phenomenon that the critical state line (CSL) for soils is parallel to the normal compression line [5,6], at least in conventional critical state soil mechanics. In this work, DEM simulations are used to establish a critical state line, and after ascertaining whether it is parallel with the NCL, to clarify what determines the slope and position of the CSL. In additional, we examine the micro mechanics behind the CSL, in particular the co-ordination number and the importance of the 'smallest particles', which must be rigorously defined.

DEM model
The simulation results presented here were performed using the software PFC [7]. For computational efficiency all particles are modelled using spheres, and gravity is neglected. All simulations are performed using a cylindrical triaxial sample, enclosed vertically by rigid platens and laterally by a flexible, faceted cylindrical membrane, shown in Fig. 1. The membrane consists of 4320 facets and the vertices can move independently from one another. Essential simulation details are provided in Table 1. The initial sample (before isotropic compression) consists of approximately 700 randomly packed particles of 2 mm diameter; although the sample is compressed and crushed (generating greater quantities of smaller particles) before any shear tests are performed.
The work uses a simple crushing model and realistic particle strengths, which the authors have previously used to produce normal compression lines [3,8,9], which exhibit the same slope as observed experimentally for the sand. This crushing model uses the average octahedral shear stress as the particle fracture criterion: which is calculated from the average principal stresses in a particle (σ 1,2,3 ), in turn calculated from contacts with neighbouring bodies. The strength data used here (and previously) is for a silica sand. From particle crushing tests [10], it was found that the size-effect on strength, i.e. the value of b in Eq. (2) was approximately 1 [10]. For a given size of particle, the strengths are distributed according to a Weibull distribution, defined by a characteristic value q 0 and a modulus m. The strength q 0 is a value such that 37% of particles of that size are stronger, and is numerically similar (and proportional) to the mean value. The modulus m defines the variability. For the silica sand, the modulus was found to be 3.3 [10]. The size-effect on strength for the simulations can therefore be written as: where for example, q 0 = 37.5 MPa for d = 2 mm. The particle breakage model replaces broken particles with 2 smaller fragments whilst obeying conservation of mass. The new fragments are placed within the boundary of the broken particle, aligned in the direction of the minor principal stress and therefore overlap. The fragments therefore move apart immediately following breakage (which requires a number of timesteps to be completed to allow energy dissipation, via friction and mechanical damping), details of which are given in [3], along with comparisons of alternative breakage models.
This initial sample is first isotropically compressed to establish a compression line, and shear tests are performed from various points along the loading path. The sample is also unloaded from different stresses, and further shear tests performed from along overconsolidated states.
During normal compression, the simulation proceeds by applying isotropic stress increments (100 kPa) using the wall servo-controls. Once a stress is applied, all particles are checked and allowed to break if their strength is exceeded. At this point a number of timesteps are completed to allow the overlapping fragments to move apart. If this breakage results in a drop in the applied isotropic pressure, then the stress increment is reapplied. Once a stress is sustainable without any further breakages, the simulation progresses to the following increment.
The shear tests are strain-controlled. The upper platen is gradually accelerated (to 0.01 m/s) then decelerated, applying an axial strain increment of 0.1%. During this process the membrane servo-control acts to adjust the radial stress: the confining pressure is adjusted to ensure either constant σ 3 , p′ or volume, depending on the stress path. Once the strain increment has been applied, particles are allowed to break, and the radial pressure is then reapplied if necessary. Once this process is complete, the subsequent strain increment is applied.
To enable the extensive number of simulations required in an acceptable timeframe, the initial sample contained only 723 particles (Fig. 1a). This is due to the crushing that occurs in nearly all shearing simulations, resulting in a wide range of sizes and much greater quantities of particles (Fig. 1c). A feature of the simulations is that there is no comminution limit-i.e. no lower limit to particle size. This allows the unadulterated evolution of a (fractal) particle size distribution. A consequence of this however is that the numerical timestep used in the simulations (influenced by d min ) becomes very small. Thus, to enable the number of simulations required whilst not imposing artificial conditions, the initial number of particles is small, which helps minimise the calculation time once extensive crushing has occurred.

Normal compression-and critical state lines
Isotropic compression results and two unloading lines from 20 and 30 MPa are presented in Fig. 2 on log e-log p′ axes, also showing the points from which the shear tests were performed. From most of these points, 3 types of triaxial test were performed: constant-σ 3 , constant-p′, constant-volume. However, it was not possible to perform constantvolume tests at the lowest p′, due to dry liquefaction (i.e. without pore fluid); and it was not computationally feasible to perform all of the constant-σ 3 and constant-p′ tests at the highest stresses due to extensive crushing caused by very large magnitudes of q and p′.
The critical states obtained from these tests are shown in Fig. 3. Triaxial simulations were run until reaching axial strains of between 25 and 35%, typically 30%, at which point all of the tests exhibited no (or in a few cases, negligible) volume or stress change. The ultimate states  from all tests were taken as critical states. Additional details of these tests, including select stress-strain results and stress-dilatancy plots are available in [11]. Fig. 3(a) includes the phase transformation points [12] which are often found to coincide with the CSL [5]. (The phase transformation points are the points at which in undrained tests, the behaviour changes from contractile to dilatant, and the rate of change of p′ reverses.) The critical states reveal a clear CSL with M ≈ 0.7. Fig. 3(b) shows the critical states on log e-log p′ axes, and at high stresses (in this case, approximately > 10 MPa) there is a unique CSL, parallel to the NCL. At lower stresses (< 10 MPa), the CSL appears nonunique, with the overconsolidated tests reaching a lower CSL, in agreement with experimental data [13].
The parallel nature of the CSL and NCL at high stresses suggest that the principles that govern the slope of the normal compression line also apply to the critical state line. The b value for the silica sand, as given in Eq. (4), is 1. Using this b in Eq. (1), predicts a slope of −0.5 (using p′ as the stress variable). Both lines demonstrate agreement with this in Fig. 3 Fig. 4(a) shows progressive PSDs obtained from an example triaxial shear simulation, which shows how the grading changes during shearing and that the rate of change reduces noticeably at large strains. Fig. 4(b) gives the ultimate PSDs obtained from the end of all triaxial shear tests. This conventional plot indicates the effects of stress level, with example values of critical state mean stress annotated. It is clear that shearing under increasing stresses leads to more developed grading curves. The same data is replotted in (c) in terms of number on logarithmic axes. For a discrete PSD (i.e. without continuous sizes), a fractal distribution will appear linear according to:

Fractal particle size distributions
where D is the fractal dimension [9,14]. It is evident from Fig. 4(c) that for critical states under increasing stress, the PSDs evolve towards a fractal distribution with a fractal dimension of around 2.5 (indicated by the dashed line), just as occurs during normal compression [14]. In theory, a fractal distribution is infinite, however in reality a fractal PSD must be limited [15], bounded by an upper and lower particle size. The largest particle size in all cases here is 2 mm, whilst the smallest fractal size depends on the extent of crushing, whereby more crushing leads to a larger range of fractal sizes. Hence in Fig. 4(c), the smallest size in each simulation is still emerging (a gradual process), and does not yet adhere to the ideal fractal shown by the dashed trendline.

Co-ordination number
One interesting observation from the simulations was that both the normal compression line and critical state line appear to possess unique coordination numbers at high stresses. The average coordination number is plotted as a function of p′ in Fig. 5 for both normal compression and critical states. During compression but before yield (approximately elastic behaviour), there is a correlation between coordination number and stress, the coordination number then drops with the onset of crushing. Beyond yield, the average coordination number during compression then remains approximately constant at around 5.6. A similar trend can be observed for critical states-large variation at the lowest stresses, before crushing becomes prominent after which the average coordination remains approximately 5.1 for states on the high-stress CSL, independent of stress level (and therefore PSD). This suggests that the average coordination number on any plastic loading line in log e-log p′ space is related to the stress ratio η (=q/p). What this also shows is that for any shear test, shearing from any point on the NCL causes the average coordination number to reduce from 5.6 to 5.1.

Micro mechanics of critical states
The kinematics behind the compression law in Eq. (1) are detailed in McDowell and de Bono [3]. The key assumptions were quantitatively examined in [9], and the theory shown to be valid for a range of real sands [4]. It will now be investigated whether the law is applicable to the critical state line. One of the assumptions was that a fractal PSD emerges during compression, which also appears the case for critical states under increasing stresses. The remaining assumptions will now be summarised, whilst full details and the kinematics behind the law can be found in [2,3,9]. In deriving the compression law in Eq. (1), it was assumed that the strength of the smallest (strongest) particles is proportional to the applied stress: which accounts for the plastic hardening of the soil, as fragments become stronger with decreasing size. If a fractal distribution evolves through the fracturing of the smallest grains, and those smallest grains are in self-similar configurations as the smallest particles become stronger and stress increases, then the strength of the smallest particles loaded in those configurations must be proportional to the applied macroscopic stress level. Recalling the particle size-hardening law in Eq. (4) therefore provides a link between applied stress and the size of the smallest particles: For a discrete particle size distribution, using the smallest particle size d s in Eq. (5) gives an expression for the number of the smallest particles: and the collective volume of the smallest particles can therefore be expressed: The fundamental assumption invoked in the derivation of the compression law given in [2] is that the volume of the smallest particles is directly proportional to the volume of voids, and therefore the voids ratio: The rationale for this is that as the smallest particles become smaller and fill the available voids, the available voids become smaller and the typical void space must be more closely related to the volume of one of those particles, as opposed to, say, the largest particle volume or the median one.
Substituting (7) into (9), together with the above assumption (10), is therefore what gives the slope of the compression line, which can be written in the form: ∝ − e p 0.5 (11) Considering that the CSL is parallel with the NCL, with a constant critical stress ratio, it seems logical that the same relationship between the volume of smallest particles and voids ratio (10) also exists on the critical state line. This will now be assessed, i.e. whether the voids ratio is proportional to the collective volume of all the smallest particles on the critical state line.
The theoretical derivation by McDowell [2] and substantiated by McDowell and de Bono [3] would work for any constant stress ratio. Thus, there should exist in log e-log p′ space, a set of parallel lines, each corresponding to a unique stress ratio, and the rate of separation of these lines as stress ratio increases governs the shape of the yield surface in stress space.
However, as shown in Fig. 4 the 'absolute' smallest particle size is not necessarily the smallest fractal particle size. (Soil does not obey a perfectly discrete hierarchical splitting model, with all particles of one size breaking at one stress; it is a gradual and continuous process.) So what counts as the smallest particles in Eq. (10) is non-trivial for a gradually evolving and expanding PSD. In general, there are only ever a few (< 10) particles of the smallest size, in some cases only 2 particles-which would clearly be inappropriate [9]. Any definition of the 'smallest' particles must contain information about the overall pore space and its interconnectivity.
It was found previously [9] that a suitable definition of the smallest particles are those which have minimal contacts-which should be a characteristic of any 'locally' smallest particles. For a particle to have a very high coordination number, it must be surrounded by smaller particles, and therefore is clearly not the smallest.
Considering random packing, the minimum contacts that a particle may have whilst mechanically loaded is 4, and this was used as the criterion to define the 'smallest' particles when investigating the above assumptions during normal compression [9]. However, it is of course possible for a small particle to have more contacts-the theoretical maximum number of contacts a spherical particle may have with equalor-larger particles is 12 (hexagonal close packing); higher contact numbers are only possible with smaller neighbouring spheres. For dense random packing however (such is the case here), this number will be much lower.
In the present analysis, it is found that using 5 contacts to define the smallest particles (and calculating their cumulative volume) leads to the best agreement in Eq. (10). It is anticipated that this difference is related to the boundary conditions-the previous analysis [9] used a rigid oedometer, with no lateral strain, and a non-uniform distribution of radial stress, while the current work uses a triaxial sample with a flexible membrane and uniform distribution of radial pressure. (Nonetheless, using a criterion of 5 contacts in the previous analysis still demonstrates acceptable agreement.) The term 'smallest particles' will now be used to refer to all particles in a simulation that have 5-or-fewer contacts. The 'volume of the smallest particles' refers to the collective volume of all such particles with 5or-fewer contacts. Fig. 6(a) shows the voids ratio as a function of the volume of these smallest particles, denoted V sm . The volume V sm used here is the volume of smallest particles normalised by the total volume of solids (constant throughout all tests), to give the percentage by volume of these smallest particles (avoiding unsightly small values). A direct proportionality exists between e and V sm , in accordance with and reaffirming the assumption in the compression law, and showing that the same relation exists at critical states. This therefore appears to explain why the CSL is parallel with the NCL, and demonstrates that the same principles apply to CSL-whereby the slope is a function of the sizeeffect on particle strengths. The voids ratio at any critical state is determined by the extent of the fractal PSD (i.e. the range of particle sizes), which itself results from crushing caused by the applied stresses. Fig. 6 shows e ∝ V sm applies to all states that reach either the plastic normal compression or critical state lines, and this includes some samples sheared from 'virgin' states, all normally consolidated states, and (lightly-overconsolidated) samples that are unloaded to lower stresses and then sheared. So samples with completely different initial PSDs reach the same e ∝ V sm line through shearing and the resulting particle crushing. For example, the compacted and over-consolidated samples sheared from, say 15 MPa, which have initial different voids ratios and PSDs (and the latter on an unloading/reloading line), reach the same critical state, with almost identical e and V sm .
An interesting observation in Fig. 6(a) is that the two sets of soil states have different coefficients of proportionality. It might have been anticipated that, if the assumption e ∝ V sm is true, then any change in voids ratio resulting from triaxial shearing would simply be due to this volume of smallest particles changing. It was shown earlier that the average coordination number during normal compression and at critical states both appear to be independent of stress level and PSD, implying that on the CSL (lower coordination number), particles are packed less efficiently together, which would explain the different coefficients of proportionalities observed in Fig. 6(a). It is now worth exploring this phenomenon further.
The equation e ∝ V sm allows the constant of proportionality to be different for each constant stress ratio. For each macroscopic stress ratio, the contact distributions will be different, giving different average co-ordination numbers, and the plastic strain increment vectors will be different at different constant stress ratios on the state boundary surface when crushing occurs. Fig. 6(a) implies that for each constant stress ratio in between the critical state stress ratio and the isotropic condition, there exists a line for each stress ratio. If there exists a unique e ∝ V sm line for each stress ratio, as well as a unique compression line (e ∝ p -0.5 ) for each stress ratio in log e-log p′ space, then the only way of satisfying both conditions would be if a unique relationship existed between V sm and p′ for all plastic states on the loose side of critical. This concept is tested in Fig. 6(b) and is shown to be absolutely the case. What this, in effect means is that Fig. 6(b) represents all plastic states on the loose side of critical and is therefore a two-dimensional representation and micro mechanical interpretation of the state boundary surface. Fig. 6(a) only shows states on the linear, high-stress NCL and CSL. Fig. 7(a) provides a fuller picture by including additional states: compacted states before yielding, OC states on the two unloading-lines, and states on the low-stress critical state lines. Fig. 7(b) provides a schematic to illustrate the behaviour. Upon normal compression the state moves along a virgin-loading line, where the voids ratio decreases only slightly with the volume of the smallest particles, with no proportionality observed between these two quantities. Upon reaching the yield stress there then appears a direct link between V sm and e, and both decrease in proportion with each other as the macroscopic stress increases. Upon unloading, the state moves back along an unloading-line where the voids ratio changes only slightly (shown by a dashed line), similar to the initial virgin-loading line.
A similar relationship exists for critical states-there also is a direct proportionality between V sm and e, however this lies to the right of the normal compression states. The relationship e ∝ V sm only applies to states on the high-stress CSL (parallel to the NCL in log e-log p′ space). Critical states at low stresses do not conform to the relationship e ∝ V sm . These states form a separate, flatter line, above the corresponding virgin/unloading-line. Fig. 7(b) shows a remarkable mirror image to the normal compression and critical state lines in log e-log p′ space.
Two constant-p′ simulations will now be examined more closely-one dilatant (compacted to 8 MPa), and one contractile (normally consolidated to 20 MPa). The paths of these simulations are shown in Fig. 8(a), in terms of e and V sm . Both simulations start on the (black) normal compression line and end on the (red 1 ) critical state line. This behaviour is generalised in Fig. 8(b). The compacted sample, which initially does not lie on the directly-proportional e ∝ V sm line, but instead on the virgin loading line, increases in volume and moves upwards to reach the 'compacted' critical state line. The NC sample, at high stresses on the NCL, decreases in voids ratio, moving downwards until reaching the critical state line.
During the 8 MPa test, in which there is negligible crushing, e increases as the sample dilates. Fig. 8(a) shows that V sm increases slightly, which might be expected from volumetric dilation-as the sample expands, particles become less-closely packed, and in general will therefore have fewer contacts. This means that more particles will be considered as the smallest, so V sm increases. This is reflected in Fig. 5(b) which showed that the coordination is lower on the critical states line.
During the 20 MPa constant-p' test, there is significant crushing and contraction, and as shown in Fig. 8(a) V sm does not decrease significantly. Ordinarily, particle crushing accompanied by a decrease in voids ratio would lead to V sm decreasing-just as occurs during normal compression (Fig. 7). The fact that V sm remains constant despite crushing and contraction suggests that there is a counteracting effect which contributes to V sm . Considering that the coordination number reduces during shearing (together with observations from the 8 MPa dilatant simulation), it would appear that particles are packed less-efficiently at the critical state. In other words, although smaller particles are produced by crushing, which enable the sample to occupy a smaller volume, the particles are relatively looser, which is supported by the observation that the average coordination number reduces from 5.6 to 5.1.
The observation that V sm remains constant during the 20 MPa simulation is consistent with the unique V sm ∝ p -0.5 relationship shown earlier in Fig. 6(b). This unique 'plastic state line' is replotted schematically in Fig. 9. According to this figure, for a normally-consolidated contractile sample, there should be no change in V sm as the state moves from the plastic NCL to the CSL while p' stays constant. The constant-p′ simulation at 8 MPa however is not initially on the plastic NCL, and so the initial state is not located on the unique (V sm ∝ p -0.5 ) plastic state line. This is also illustrated in Fig. 9, where elastic loading states are located below the unique V sm -p line. As shown, when such dilatant Fig. 6. Voids ratio as a function of the volume of smallest particles V sm (a), and V sm as a function of mean stress. 1 For interpretation of color in Fig. 8, the reader is referred to the web version of this article.
samples are sheared at constant-p′, the state moves upwards as V sm increases, towards the unique plastic state line.
This diagram can also explain the observed behaviour of constantvolume tests. For contractile constant-volume tests, such as that at p′ = 25 MPa, the mean stress reduces, so according to this figure the state moves leftwards up the plastic state line, and so V sm increases. This behaviour can be visualised in Fig. 8 as a path at constant e directly right from the NCL to the CSL. Just how V sm increases is due to particle breakage occurring at a constant volume-shearing induces stresses great enough to cause significant particle crushing, which causes p′ to reduce to maintain constant volume. However, the fact that volume change is prevented indicates that the sample becomes relatively looser, analogous to e min decreasing, and means that in general particles will have fewer contacts and therefore more will be classed as the smallest, causing V sm to increase. This again is reflected in a reduction of the average coordination number, which reduces from 5.6 to around 5.1.
Equally, for dilatant constant-volume tests, in which negligible breakage occurs, the mean stress increases to suppress dilation. As shown in Fig. 9, this is represented by a path that moves rightwards and down the V sm ∝ p -0.5 line, and a reduction in V sm . An example of this behaviour is the OC test at 5 MPa: p′ increases to over 13 MPa; and with only limited crushing, this means that the particles are more tightly packed together and have more contacts. This means fewer particles have ≤5 contacts and are considered as the smallest, and therefore V sm  decreases. Fig. 9 is also consistent with the observed behaviour for constant-σ 3 triaxial tests (indicated in the figure), as well as normal compression-in which the state simply moves down the plastic state line.

Conclusions
A cylindrical sample of spheres was isotropically compressed to establish a normal compression line, and then a large set of shear tests were performed from various states to establish the critical state line, which was parallel with the NCL at high stresses, in accordance with general experimental behaviour. From analysis of the particle size distributions at critical state, it was shown that the crushing resulting from shearing contributed to the emergence of fractal PSDs with a fractal dimension of around 2.5, similar to what occurs during normal compression. This observsation, coupled with the parallel nature of the CSL and NCL, suggested that the same principles behind the compression law could be used to reveal the micromechanics behind critical states.
An assumption of the compression law is that on the NCL, the available void space is proportional to the cumulative volume of the smallest particles. What constitutes the 'smallest' particles has previously been examined, and the best measure was found to be using the coordination number to define the smallest particles. In the simulations presented here if the smallest particles are defined as those with 5 or fewer contacts, then the total volume of such particles appears to be directly proportional to the voids ratio for both isotropic compression and critical states. Thus, the parallel critical state line is a result of the fractal crushing of particles, where the slope is a function of the sizeeffect on particle strength. The voids ratio at a critical state is determined by the extent of the fractal PSD, caused by crushing due to the applied stress.
The value of 5 is, however, unimportant. What is important is the dynamic definition of the smallest particle size in a constantly evolving PSD. Using a minimal value of coordination number seems a logical way of defining the (locally) smallest particles in the sample, and enables an explanation from a micro mechanic point of view of the slopes of the NCL and CSL as well as the observed behaviour of compacted, normally consolidated and overconsolidated sands. This is a scientific approach to defining the smallest particles, which includes information about the surrounding pore space; more primitive definitions might use the physically smallest particle(s), which may have no load-bearing ability or there may only be 1 such particle, as soil does not crush hierarchically at discrete stresses.
The proportionality between current voids ratio and the volume of the smallest particles was valid for all states on the plastic, high-stress normal compression or critical state lines, but the coefficient of proportionality is different for each stress ratio. It was shown however that a single, unique relation exists between the volume of the smallest particles and the current applied stress, for any stress ratio. When plotted graphically, this single relation represents all states on any plastic loading line (e.g. CSL or NCL), and therefore could be used to explain the behaviour of any sample subjected to triaxial shearing. This newly established unique relationship between V sm and p′ is effectively a two-dimensional interpretation of the state boundary surface and future work will examine the evolution of V sm along various stress paths at constant stress ratio.