Improving constant-volume simulations of undrained behaviour in DEM

In order to simulate undrained conditions using the discrete element method, a constant sample volume is often assumed. There are well-recognised problems with these constant-volume triaxial simulations, particularly of dense samples, which inhibit quantitative comparison with laboratory experiments. In this paper, four possible explanations for these problems with conventional constant-volume simulations of ideal spherical particles are explored, each of which has a physical basis: particle crushing, the presence of highly compressible air within the sample, or the reduction in stiffness due to particle surface asperities or non-spherical particle shapes. These options are explored independently and in combination through implementation in the open-source LAMMPS code. In situations where a significant amount of particle crushing occurs, it is important to incorporate this in the simulations so that stresses are not over-estimated. There is experimental evidence that irregular particles have lower Young’s moduli than the Hertzian spheres often used in DEM. In the absence of particle crushing, the most effective method to achieve more realistic stress–strain responses is to reduce the particle shear modulus substantially. This approach has the added computational benefit of enabling an increase in the simulation time-step.


Introduction
Soil is a complex multi-phase material consisting of solid, liquid and gas. Undrained tests permit the behaviour of soil to be investigated from which the pore fluid does not have sufficient time to escape when subjected to load. Excess pore pressure is generated during shearing under undrained conditions. This excess pore pressure controls important soil responses such as liquefaction: the complete loss of soil strength and stiffness [6,25]. Many researchers, e.g. [3,4,31], have studied undrained soil behaviour using laboratory testing. These tests enable us to understand the macro-scale responses of soil such as the stress-strain behaviour but cannot give any information on the dynamic changes occurring at the micro-scale that cause the observed macro-scale response. In recent years, the discrete element method (DEM) [14] has become very popular in geomechanics research due to its ability to capture the macro-scale response of soil while enabling investigation at the micro-scale [34]. Generally, there are two modelling approaches used by researchers to simulate undrained tests using DEM. The first approach involves coupling the DEM code with a suitable fluid-solving code, often computational fluid dynamics [39,50]. However, this adds complications to the simulation and increases the computational cost.
In the alternative 'constant-volume' approach, the sample volume is maintained constant throughout shearing and the excess pore water pressure is estimated as [49]: 3;o is the initial confining effective stress at the start of shearing and r 0 3 is the minor principal effective stress at every subsequent time-step in the DEM simulation. Constant-volume (CV) simulations assume that the soil is fully saturated and water is incompressible. This approach has been widely adopted, e.g. [16,49]. The constant-volume method has the advantage of computational simplicity. However, some problems arise when shearing dense samples, stemming from the generation of unrealistically high stresses [19]. Indeed, stresses [ 10 MPa are often generated in CV simulations which invalidate the underlying assumption of the incompressibility of water, e.g. the volume of water compressed to 10 MPa is reduced by 0.5%, given a bulk modulus of 2.2 GPa [19]. The assumption of point contact in DEM [37] is also violated, e.g. overlaps reached almost 20% (normalised by mean diameter) by the end of the simulations presented by Hanley et al. [19]. Other challenges include strain rate sensitivity during simulations [19], even at very low strain rates, unrealistically high small-strain stiffness, and the failure of the simulations to capture the volume changes observed during undrained laboratory shearing [7,33,41].
These problems mean that an alternative to the constantvolume method should be sought which retains the method's computational efficiency but without the unphysicality for dense soils. The aim of this paper is to establish a computationally efficient and physically justifiable alternative to the conventional constant-volume method using perfect spheres in DEM. Four alternatives are hypothesised to address the shortcomings of the constant-volume method described above: 1. Particle crushing causes a reduction of the stresses upon shearing. 2. Air greatly increases compressibility of the pore fluid, causing some changes in sample volume. 3. Particle surface asperities reduce the initial contact stiffness. 4. Adopting Hertzian spheres to represent non-spherical particles gives an overly stiff response. This can be corrected by reducing the particle shear modulus by a factor obtainable from uniaxial compression of single particles.
These options are explored independently and in combination through implementation in the open-source LAMMPS code [36]. Based on this study, recommendations are made to improve quantitative agreement with laboratory data.
2 Theoretical background 2.1 Influence of particle crushing Particle crushing often occurs during shearing or compression of real sands [1,30]. DEM simulations of triaxial shearing have shown that particle crushing causes the peak stresses to reduce, the volumetric response to become more contractive and the position of the critical state line to shift in e-log(p 0 ) space [5,15,20]. Even though it can have significant effects, particle crushing is often disregarded in DEM simulations for two reasons: (1) the necessity to simulate fine particles post-crushing is very computationally expensive; (2) the assumptions that are made for reasons of computational tractability can be unphysical, e.g. accepting the loss of solid volume from the simulation, imposing a comminution limit or predefining highly idealised fragment size distributions. A recently developed crushing model from the literature was adopted for this study [20]; the reader is directed to that paper for a detailed discussion of the development and implementation of the model. In summary, a particle is deemed to fail when any contact force acting on the particle exceeds a predefined crushing force. These crushing forces are experimentally measured from uniaxial compression of single particles [32]. The inherent variability in these crushing forces is captured using Weibull statistics in the model. Upon failure, the particle's radius is reduced so that contact is lost with all surrounding particles and its crushing force is increased. Fine particles are inserted into the void space to conserve solid volume. Particles can no longer fail once a comminution limit has been reached. This particle replacement approach to particle crushing has two major advantages compared to the agglomerate breakage approach: the sizes of daughter particles are not predefined and the fine fragments are not simulated from the start of the simulation which is computationally efficient [20].

Influence of air: bulk modulus of water-air mixtures
The pore fluid is not directly simulated in this research. Instead, the compressibility of the pore fluid is captured by allowing the total volume of the periodic cell to vary based on the theory presented in this section. Even though the pore fluid is not directly simulated, some of its properties can be estimated by subtraction, e.g. its volume must equal the total cell volume minus the solid particle volume. One cause of unrealistically high stresses may be the assumption of perfect saturation; in physical experiments, samples contain a small fraction of air. Before shearing begins, a B-test is usually performed to approximate the degree of saturation of the soil. When B = 1 or 100%, the soil is fully saturated. Typically, B = 90-95% during a physical undrained triaxial test, meaning that 5-10% of air is present in the soil sample even though it may still be considered fully saturated. Since air is highly compressible, the presence of a small percentage of air is likely to be influential.
The increment of pore pressure Du is the difference between the increments of total stress Dp and mean effective stress Dp 0 : The mean effective stress and deviator stress can be expressed in terms of principal effective stresses. Since the volume change is uniform in the radial direction, i.e. Dr 0 2 ¼ Dr 0 3 : For an undrained triaxial compression test where the intermediate stress ratio b = 0, the loading path is equal to Dp ¼ Dq = 3 . Therefore, This change in the minor principal effective stress can be computed from the interparticle contact force data in DEM. Equation 5 is the definition of the bulk modulus, assuming that the soil particles are incompressible: Vg, the product of porosity and sample volume, is the volume occupied by pore fluid, K f is the bulk modulus of the pore fluid and DV is the volume change of the soil sample, i.e. the volume change of the pore fluid. Reorganising Eq. 5, we get: V current is the total sample volume updated at every timestep during the simulation and V particles is the fixed volume of the solid particles within the soil sample. The minus sign in Eq. 6 indicates that the sample volume is allowed to expand when Du is negative (or Dr 0 3 is positive) during shearing. Conversely, a positive Du in Eq. 6 leads to sample contraction, capturing compression of the entrained air that takes place in the physical test.
Knowing the bulk modulus of the pore fluid, K f , the volume change at every time-step may be calculated from Eq. 6. K f is given by [44,46]: K w is the bulk modulus of water, S r is the degree of saturation of the soil sample and P a is the absolute fluid pressure: the sum of atmospheric pressure and excess pore pressure. The assumption being made is that the degree of saturation is high, the pore fluid is homogenous and the air exists in the pore water in the form of well-distributed bubbles [46]. The dissolution of air into water at high pressure has been neglected. Figure 1 shows the variation of K f with P a according to Eq. 7, taking K w as 2.2 GPa. When the degree of saturation is 100%, the bulk modulus of the pore fluid is equal to the bulk modulus of water. K f is substantially reduced by the presence of a small percentage of air at low to moderate pressures relevant to laboratory soil testing.

Influence of particle surface asperities:
rough-surface contact model Another explanation for the unrealistically high computed stresses is the perfectly smooth nature of the interparticle contact. The presence of surface asperities on real particles reduces the contact stiffness compared to smooth spheres during the initial phase of loading. The influence of surface asperities at interparticle contacts is discussed by [18,47,48]. In this research, a DEM contact model developed by Otsubo et al. [35] was adopted which includes crushing of asperities. This contact model was created as an extension of a previous model developed based on single-particle compression tests which includes surface roughness, S q , and hardness [11]. Otsubo et al.'s model [35] includes three regimes: asperities dominating, a transitional regime and Hertzian contact. d T1 and d T2 are the threshold contact displacements at N ¼ N T1 and N T2 , respectively, as shown in Fig. 2. At interparticle contact overlaps less than d T1 , the contact response is dominated by crushing of surface asperities. This is controlled by two

Influence of using Hertzian spheres to represent non-spherical particles
The presence of asperities on the surface of the particle is expected to reduce the initial contact stiffness. However, experimental data for uniaxial compression of individual particles [8,10] suggest that the overall particle shape may have a much more significant effect on the load-deformation response. For an irregular silica gravel particle, the experimentally measured Young's modulus was found to be 20 times less than the Young's modulus of an equivalent sphere calculated from Hertzian mechanics [9]. This experimental result must be reconciled with other experiments [45] which show that increasing the particles' angularity while keeping all other inputs constant gives a stronger bulk stress-strain response. The simple addition of rotational resistance, to include some degree of shape irregularity, has a similar effect in DEM simulations [24]. The reason for this apparent inconsistency is that nonspherical particles do not behave according to Hertzian mechanics of spheres, the adoption of which gives an overly stiff response. In the absence of a better understanding of contact mechanics for non-spherical particle shapes, the contact forces and stiffness calculated using Hertzian mechanics of spheres can be corrected by reducing the Young's modulus of the particle (which is, in reality, an irregular sand grain). This reduction could be calibrated to obtain the correct stiffness for an individual particle under uniaxial compression. The reduction in the stiffness of the bulk sample is purely a result of reducing the stiffness of each individual particle comprising the sample. The results in [9] imply that the bulk stiffness of an assembly of irregular particles would be substantially lower than for perfect Hertzian spheres composed of the same material.

Code implementation
All of the simulations were run using a version of opensource, MPI-parallelised LAMMPS code [36]. The approaches by which the four hypotheses described in Sects. 2.1-2.4 were implemented in LAMMPS are discussed in Sects. 3.1-3.4, respectively.
3.1 Influence of particle crushing The particle crushing model described in Sect. 2.1 was developed by Hanley et al. [20]. The implementation of this model in LAMMPS was carried out as part of [20]; the reader is referred to that paper for details.  Fig. 2 Schematic of the rough-surface contact model used in this study, based on [35] 3.2 Influence of air: bulk modulus of water-air mixtures The flowchart in Fig. 3 shows the implementation of the equations described in Sect. 2.2 in LAMMPS. In summary, this version of LAMMPS contains a servo-control algorithm for periodically bounded samples which requires the computation of principal effective stresses from interparticle contact forces. Dr 0 3 is found as the change in minor principal effective stress between successive time-steps. This is equal to the negative change in the excess pore water pressure, ÀDu, i.e. Equation 4. Du is accumulated and is added to atmospheric pressure, P atm , to find P a , the absolute pore fluid pressure. Two limits were imposed in the implementation of the equations described in Sect. 2.2. The first limit was imposed on P a to avoid the attainment of low pressures which would cause fluid vaporisation in reality, and beyond u += ∆u Stored ∆ε V *= Dissipation factor Stored ∆ε V += ∆ε Vlim Stored ∆ε V -= ∆ε Vlim ∆σ 3 ' = -∆u P a = u + P atm Stored ∆u = P lim -P atm Revise increment of ∆u P a = P lim P a < P lim Find K f using P a Stored ∆ε V += 0.5 * (∆V/V current ) |Stored ∆ε V | > ∆ε Vlim that, the attainment of non-physical negative pressures. A limit of P a ! 0:25P atm was imposed. This predefined limit is the lowest permissible absolute pressure. The upper bound on P a is the confining pressure, r 3 , at which the effective stress terms become zero. P a is used to calculate the bulk modulus of the pore fluid, K f , using Eq. 7. The volume change of the sample on that time-step is given by Eq. 6, assuming a discontinuous jump from time-step t to the time-step t þ Dt. If the time-step Dt were divided into infinitely many smaller time increments, the required volume change would be half of that given by Eq. 6. This halved volume increment corresponding to a continuoustime scenario was implemented in the code. The derivation is given in ''Appendix''.
The second limit was imposed on absolute volumetric strain on each time-step of 2.3 9 10 -9 . This value was determined through trial and error so that the absolute volumetric strain during shearing did not exceed the volume of air present in the pore fluid (5% of the pore fluid). Without this condition, huge values of DV could be achieved in single time-step when very stiff particles are used and hence Dr 0 3 could be unreasonably large. This would lead to instability of the proportional servo-controller. Volumetric strains beyond the limit are stored until the following time-step. The stored volumetric strain term is multiplied by a dissipation factor of 0.9999 to prevent the accrual of huge stored strains and ensure that critical state is eventually reached. Using K f and Du, the volume change is calculated at every time-step using Eq. 6.

Influence of particle surface asperities:
rough-surface contact model The rough-surface contact model described in Sect. 2.3 was developed by Otsubo et al. [35]. The implementation of the model in LAMMPS was carried out as part of [35]; the reader is referred to that paper for details.

Influence of using Hertzian spheres to represent non-spherical particles
The Young's modulus of an irregular silica gravel particle measured experimentally by Cavarretta and O'Sullivan [9] was 20 times lower than that of an equivalent Hertzian sphere. In this research using Hertzian spheres, the shear modulus of the individual particles, G, was reduced by a factor of 20 to compensate for this known disparity between experimental data and the Hertzian predictions for spheres.
An added benefit of this approach is an increase in the stable simulation time-step. For the nonlinear Hertzian contact model, the critical time-step based on the Rayleigh wave velocity for the system has the relationship Dt c / G À0:5 [43]. Thus, reducing G by a factor of 20 increases the time-step by a factor of ffiffiffiffiffi 20 p % 4:5, substantially reducing a simulation's run-time without compromising its stability.

DEM model preparation and simulation plan
Cubic granular specimens of 10 9 10 9 10 mm 3 were created which contained 28,309 spheres. The grading used was not representative of any specific sand; particle diameters varied between 0.1 and 1 mm ( 50 = 0.516 mm), with = 3.004 and = 0.573 indicating a poorly graded sand. To eliminate boundary effects [22] and to ensure homogenous deformation [13], periodic boundary conditions were adopted for this study. The particles were placed randomly within these periodic cells without initial overlaps using a MATLAB code. Parametric studies confirmed that these periodically bounded samples contain sufficient particles to ensure a representative element volume [22,27]. The initial particle positions were imported to LAMMPS before each sample was isotropically compressed by moving the boundaries under stress control to achieve a specified confining pressure of 150 kPa (Fig. 4).
The friction coefficient was set to zero during the sample preparation process to generate dense samples. This friction coefficient was increased to 0.25, based on [23], before shearing each sample at a fixed strain rate of 1 s -1 . This strain rate ensured an inertia number less than 3.6 9 10 -5 throughout shearing: lower than the limiting value of 7.9 9 10 -5 proposed for quasi-static behaviour [29]. Eight triaxial shearing simulations are the principal focus of this study: four using the constant-volume method and four using the bulk modulus method. Each subset of four simulations consisted of simulations with the simplified Hertz-Mindlin (smooth) [26] and rough-surface [35] contact models, with shear moduli of 29 GPa or 1.46 GPa. These shear moduli, G, respectively represented a physically realistic value for quartz and a reduced value to capture irregularity of the particle shape based on Cavarretta and O'Sullivan [9]. The particle density and Poisson's ratio were set at 2675 kg/m 3 and 0.2, respectively. The local damping coefficient was chosen as 0.2. The bulk modulus of water, K w , was 2.2 GPa and the degree of saturation, S r , of the soil sample was 0.95. Atmospheric pressure was set at 0.1 MPa. The surface roughness, S q , was set at 0.5 9 10 -6 m for the rough-surface contact model: similar to that of an LBSA sand grain [8]. Values of d 1 ¼ 0:82S q and d 2 ¼ 1:24S q were used in these simulations, based on experimental data [18,48]. Gravity was not considered and particle crushing was not permitted in these eight simulations.
It is already well known that particle crushing reduces the peak stresses in a triaxial test, e.g. [20]. Only one simulation was run in which particle crushing was considered. This simulation used the constant-volume method, simplified Hertz-Mindlin contact model and G = 29 GPa. The parameters of the crushing model were a Weibull modulus of 4.2, a limiting comminution radius of 50 lm, and a characteristic stress (r 0 ) of 760 MPa at which 37% of the particles of characteristic diameter 1.29 mm survive. A linear trendline of the form P s d ð Þ ¼ a r r 0;d þ b, relating the probability of survival for particles of diameter d to a stress r, was assumed, based on the statistics in [32] for a quartzitic Aio sand. These parameters were identical to those in [20], except a larger r 0 value was used in this study to limit the amount of crushing that occurred. To obtain a critical state line, three supplementary triaxial drained simulations were run at confining pressures of 150, 300 and 500 kPa and two constant p 0 simulations were run at 40 and 47 MPa.

Macro-scale
Consider firstly the one simulation which includes particle crushing. This is compared to the equivalent simulation in which crushing is ignored on Fig. 5. Without crushing, the mean effective stress, p 0 , is 28.72 MPa at 15% axial strain. The inclusion of particle crushing reduced p 0 by 63% at the same strain. However, this p 0 is much higher than in a laboratory test conducted by Kuwano [28], who obtained p 0 = 1.46 MPa at 15% axial strain for Dunkerque sand with a confining pressure of 400 kPa. This experimental data are also shown in Fig. 5. The stresses generated can be calibrated using r 0 -760 MPa for this simulation-as an adjustable parameter. Using this parameter, a significant amount of crushing occurred; the number of particles increased to 52,898 (28,309 particles before shearing) when the sample was sheared to 15% axial strain. The change in particle size distribution (PSD) is shown in Fig. 6. Reducing r 0 induces more particle crushing and hence reduces p 0 . Hanley et al. [20] used r 0 = 38 MPa and obtained a huge amount of breakage at high confining pressures in drained simulations. In situations where a considerable amount of particle crushing occurs, it is important to consider this in the simulations to correctly capture the bulk behaviour.  Mean effective stress (MPa) against axial strain (%) for the constant-volume triaxial simulations with and without crushing compared with data from an undrained triaxial laboratory test for a confining pressure of 400 kPa [28] Fig. 6 Comparison of particle size distributions by number before shearing and at 15% strain for the constant-volume triaxial simulation where particle crushing is considered However, particle crushing does not fully explain the disparity between undrained laboratory tests and constantvolume DEM simulations. The amount of crushing that would be required to quantitatively match the stress-strain behaviour for a dense sample would be far more than observed in laboratory tests of sands (and would be unachievable with the type of crushing model proposed by [20]). Furthermore, undrained tests on dry spherical glass beads conducted by Cui et al. [12] observed a maximum deviator stress of 440 kPa for a confining pressure of 200 kPa and no particle crushing occurred. This indicates that particle crushing is not the only reason for high stresses in constant-volume simulations.
Considering separately the eight simulations without crushing, all simulations reached a critical state before 25% axial strain, attaining a similar stress ratio, q=p 0 , of around 0.72 as shown in Fig. 7, regardless of simulation method, contact model or particle stiffness. These plots show behaviour characteristic of dense samples: the stress ratio increased abruptly to a peak value upon initial shearing and dropped thereafter to a constant value at the critical state [33]. Figure 8 shows the variation of mean effective stress, p 0 , and deviator stress, q, with axial strain where these quantities are defined by Eq. 3. Using the constant-volume method with a Hertzian contact model and particle shear modulus of 29 GPa, the stresses generated at critical state were 95.18 MPa and 68.14 MPa for p 0 and q, respectively: similar in magnitude to the stresses reported by [19]. Changing to a rough-surface contact model allowing for asperity crushing reduced the stresses to 80.72 MPa and 58.42 MPa, respectively: a reduction of around 15%.
Switching from constant volume to the bulk modulus method was more effective, leading to a reduction of 56% (41.87 MPa and 29.89 MPa). The most effective method to achieve more physically realistic stresses was reducing the shear modulus from 29 to 1.46 GPa, which captures the effect of irregular, non-spherical particle shapes when Hertzian mechanics are adopted. This yielded a reduction of 94%, to p 0 = 5.28 MPa and q = 3.83 MPa, compared to using the shear modulus of quartz. A 20-fold reduction of G reduced the stresses to around one-twentieth of their former values. In combination, using the bulk modulus method with a rough-surface contact model and G = 1.46 GPa led to the lowest stress state among these eight simulations of p 0 = 3.58 MPa and q = 2.55 MPa at critical state. These results are quantitatively similar to experimental results [2]. The initial concavity in these stress-strain curves matches experimental and simulation data in the literature for undrained triaxial compression of dense sands [21,38].
The shear modulus of the bulk soil sample is computed as 1 3 of the slope of a plot of deviator stress against triaxial shear strain. As shown in Fig. 9, the shear modulus of the soil samples using a smooth Hertzian contact model with a particle shear modulus of 29 GPa was extremely high: 311 MPa initially at 1 9 10 -3 % strain. By considering the particle surface asperities, the shear modulus of the soil sample was considerably reduced (by 64% initially) compared to the smooth Hertzian model. The high stresses generated during shearing may be caused by this initial unrealistically stiff response of Hertzian spheres using the particle shear modulus of quartz. Using a reduced particle shear modulus of 1.46 GPa, the initial stiffness of the sample was reduced by around 85% compared to using the shear modulus of quartz particles.
The variation of volumetric strain and void ratio with axial strain during shearing are plotted on Fig. 10. The four simulations using the constant-volume method do not permit any volumetric strain throughout shearing. When the bulk modulus method is used, the samples contract slightly during initial shearing and then dilate until critical state is attained. This expansion of each sample's volume principally reflects the expansion of the air within the pore fluid as negative pore pressures develop. According to Boyle's law, the volume increases as the pressure of gas decreases within a closed system. In this implementation, it is noted that the maximum dilation achieved is smaller than indicated by Boyle's law at the limiting pressure of 0:25P atm . This is due to the volumetric strain limit imposed on each time-step during shearing, which restricted the sample dilation to 1-2% at critical state. This is similar to the volume of air within each sample: at S r ¼ 0:95, 5% of air is present in the pore fluid which corresponds to around 1.5% of the total volume of each sand sample. As stresses increase to a maximum of p 0 = 95.18 MPa for a Hertzian contact model with G = 29 GPa, dilation also increases. Figure 11a shows that the adoption of any of the proposed alternatives to the constant-volume method, i.e. the bulk modulus method, the rough-surface contact model to capture crushing of surface asperities, or the reduced particle shear modulus to correct for the use of Hertzian spheres to simulate non-spherical particles, all lead to a downward shift of the critical state line in e-log(p 0 ) space. Particle crushing at high p 0 causes a similar shift of the critical state line [20]. As the stress ratios at critical state are very similar for all simulations (Fig. 7), all points are 8°. This angle is much lower than values obtained in physical tests. The main reason for this disparity is the simulated particles are spherical and can rotate freely whereas, in a real system, interparticle movements are much more inhibited because of interlocking, both at the particle level and at the contact level where interlocking of asperities prevents free rotation.
The mean interparticle overlap normalised by mean particle radius is shown on Fig. 12 as shearing proceeds. Regardless of contact model used, the mean overlaps were below 5% in all cases. However, the largest overlap exceeded 20%. As the particle shear modulus was reduced from 29 to 1.46 GPa, the mean overlap increased slightly but remained below 5%. This is important to quantify because of the fundamental assumption of point contact in DEM [37]. The contact area, and hence the interparticle overlap, is significantly larger for the rough-surface contact model than for the smooth Hertzian contact model which agrees with Greenwood and Tripp [17]. Because the stresses are lower when the presence of air is considered (bulk modulus method), this method leads to reduced interparticle overlaps compared to the conventional constant-volume method.

Micro-scale
The coordination number, Z, is a scalar measure of fabric, i.e. internal topology, within a granular system. It is computed as in which C is the number of interparticle contacts and N is the number of particles in the sample. In all cases, Z decreases immediately upon shearing and attains a constant value of 4.2-5.0 at critical state (Fig. 13). The bulk modulus method consistently gives a substantial reduction in coordination number compared to the equivalent constant-volume simulations. This is expected as dilation is associated with a reduction in contact density within a granular assembly. The rough-surface contact model gives a slightly higher coordination number than the smooth Hertzian model. Reducing G also increases Z. Both of these methods lead to higher mean interparticle overlaps (Fig. 12) which is expected to give a small increase in  C. The mechanical coordination number, Z m , is computed similarly to Z, except those particles with zero or one interparticle contact are excluded from the calculation [42]. The trends in Z m on Fig. 13 broadly match those for Z; however, the range of Z m values at critical state of 5.4-5.7 is narrower than the range of Z values. These results quantitatively agree with Huang et al. [23].
The deviatoric fabric, u d , is the difference between the maximum and minimum eigenvalues of the second-order fabric tensor defined by Satake [40]: It is widely used to quantify the fabric anisotropy of granular assemblies. Figure 14 shows the variation of deviatoric fabric with axial strain for all simulations considered. u d is almost zero at the start of shearing of each isotropic sample. u d attains a maximum value of around 0.08 at 5% axial strain before decreasing to a stable value between 0.04 and 0.05 at critical state. For the same particle shear modulus, the bulk modulus method simulations showed slightly higher u d values at peak and critical state than the constant-volume simulations. These results also quantitatively agree with Huang et al. [29]. Reducing G decreases u d .

Conclusions
This paper was introduced with the aim of establishing a computationally efficient and physically justifiable alternative to the constant-volume method with ideal spherical particles. Four alternatives were hypothesised, each of which has a physical justification: that particle crushing substantially reduces the peak stresses upon shearing, that air greatly increases compressibility of the pore fluid, causing some changes of sample volume to occur; that particle surface asperities reduce the initial contact stiffness; and that non-spherical particle shapes reduce the sample stiffness when Hertzian mechanics for spheres are adopted for the calculation of contact forces. All of these hypotheses were explored, both independently and in combination, using triaxial compression simulations run using the LAMMPS code.
When the conventional constant-volume method was used with a Hertzian contact model and particle shear modulus of 29 GPa, p 0 exceeded 95 MPa at critical state, highlighting one of the main problems with this approach. When particle crushing was considered, the stresses were substantially reduced. In principle, the parameters controlling the degree of crushing which occurs could be calibrated to give the desired macro-scale response. However, the amount of crushing that would be required to give the correct stress-strain response would be unrealistically high for a dense sample if a constant-volume simulation with smooth spheres were chosen for the DEM. The computational expense of such a simulation would also be prohibitive.
The most effective method to achieve more realistic stresses was reducing the shear modulus by a factor of 20 (p 0 = 5.28 MPa at critical state). This captured the effect of irregular particle shape when Hertzian mechanics are adopted, based on experimental measurements of the Young's modulus of an irregular silica gravel particle [9]. This method is also computationally beneficial as the reduced shear modulus allows the simulation time-step to be increased by a factor of approximately 4.5. Furthermore, this method reduces the small-strain stiffness of the sample  to more realistic values. Adopting a rough-surface contact model (capturing the effect of asperity crushing) or switching to the bulk modulus method (enabling changes of sample volume) were both less effective and more computationally expensive than reducing G. All simulations attained similar stress ratios at critical state of around 0.72 while the deviatoric fabric was almost unaffected.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.
where Dr 0 3 3 À X 1 is the stress difference ðr 0 3;b À r 0 3;a Þ. Similarly, X 1 and X 2 therefore represent the drops in r 0 3 that result when DV 1 and DV 2 are added, respectively.
Substituting X 1 and X 2 into Eqs. 11 and 12, respectively, we get The total volume change over the time-step Dt is obtained by summation: DV 1 þ DV 2 þ DV 3 ¼ For even increments, the incremental volume change is As q ? !, it is straightforward to show that the total required volume change is equal to X DV i ! a Dr 0 3 2 This is half of the DV assuming a discontinuous jump from time-step t to time-step t þ Dt. This volume increment for a continuous time scenario was implemented in the code.