Rattlers’ involvement for possibly looser critical states under higher mean stress

The critical state of a granular material made of rounded tetrahedral particles is studied through DEM simulations of triaxial compressions. A minimum number of 7,500 particles is first obtained as the representative volume element (RVE) for the present triaxial simulations. Then, the macroscopic critical state line (CSL) is shown to be increasing in the (stress, density) space for low confining pressures. Such an unexpected behaviour is explained by the existence of a significant proportion of rattlers. Considering rattlers as voids indeed reinstates a classically decreasing CSL.


Introduction
The ability of granular materials to sustain increasing shear strains under constant volume and constant shear stress has long been recognised as one of their salient features. From a soil mechanics point of view, such a deformation mechanism has been coined as "critical state" [1] and laid the basis for countless constitutive relations [2]. Experiments [3][4][5] have shown that a critical state line (CSL) exists in the (log(p), e) plane (with p the mean pressure and e the void ratio) and that this CSL is decreasing. These observations were mostly confirmed numerically in [6][7][8], using Discrete Element Methods (DEM) with spherical particles and the Hertz-Mindlin contact model. The parametric analysis led in [6] nevertheless has shown that the slope of the CSL depends on the inter-particle friction angle, and a positive slope was actually observed therein in some cases. This striking result was interpreted considering that, for high values of inter-particle friction angle, fewer particles are needed to maintain the static equilibrium of the packing and that several particles become "rattlers", having less than two contacts. Rattlers being strangers to the force chains, they do not really form part of the solid phase, and may bias the determination of the CSL, up to getting a positive slope.
The present study pursues the investigation into a possibly increasing CSL in DEM, that would contradict experiments. The granular material adopted herein is different from [6], in terms of shape and contact model. Our numerical triaxial tests are indeed performed on a material made of tetrahedral particles interacting according to a visco-elastic contact model with friction. Note that gravity is not considered in this study. These simulations are conducted using the open source code YADE ( [9]). The material studied in this paper is constituted of identical tetrahedral particles ( Figure 1). Each particle is a clump of four spheres with the same radius R sph = 3.101 mm, being inscribed in an outer sphere of diameter D clp = 10 mm. These particles have a density of ρ = 1,111 kg m −3 so as to coincide with resin and the same holds for the other material parameters reported in Table 1. The triaxial tests are performed using six infinite walls made of acrylic which interact with the particles but not with each other. These materials were actually used experimentally by the Japanese Geotechnical Society (JGS) for an ongoing round-robin test 1 . The relevant contact parameters used in DEM were measured experimentally by the JGS (see Table 1). The internal friction angle ϕ was determined by performing sliding tests with resinous cubes against boards made of resin or acrylic. The normal restitution coefficient e n was measured by performing drop tests with resinous spheres on boards made of resin or acrylic. Lastly, the normal stiffness K n was determined on single spheres as the slope of the loading-displacement relation on the third cycle of compression tests. As for the tangential stiffness K t , it was determined by setting arbitrarily the stiffness ratio K t /K n .

Contact model
The contact model accounts for friction through the interparticle friction angle ϕ, used during the computation of the tangential force F t . It also accounts for visco-elasticity through the restitution coefficient e n , used during the computation of the normal force F n . Indeed e n can be expressed according to a damping coefficient c n [10]. This damping coefficient can be determined using the Newton-Raphson method and then be used directly in the computation of F n . Denoting u n and u t the normal and tangential relative displacements respectively, the normal and tangential forces finally read: Since the visco-elasticity already stabilizes the simulation an additional numerical damping is not necessary and was not introduced in the simulations.

Triaxial tests
The triaxial tests workflow is the following: A cloud of particles is created by selecting the positions of the clump centers randomly. In order to eliminate this source of variability, the random number generation seed was fixed so that every triaxial tests begin with the same cloud of particles.
Six walls are created around the sample, they are then moved toward each other so the sample is isotropically compacted to a nominal confining pressure P c . During this step, the inter-particle friction angle ϕ comp can be set between 0.01 • and 27 • to reach different initial packing densities. This method has already been used successfully in [11] for 2D simulations and in [12] for 3D simulations.
Once the sample is stable and the mean pressure on the walls P w verifies 0.995P c < P w < 1.005P c , the compaction is considered over and ϕ is set to the value given in Table 1. The sample stability is measured using the unbalanced force u F , defined in [9]. The chosen condition for sample stability is u F < 10 −2 . P c is then maintained on the side walls and one axial wall is moved such that a strain rate˙ ax is imposed to the  Figure 2. q and V during some S1 simulations sample. Denoting the inertial number I n = 2 × 10 −4 , low enough to ensure quasi-staticity,˙ ax is computed using the expression proposed in [13]: The triaxial test is stopped when the axial strain ax reaches 0.8. The critical state is considered to be reached at ax = 0.6, which is enough according to [3][4][5][6][7][8] and Figure 2. For any quantity s, its critical value s crit will thus be computed as its average over ax ∈ [0.6, 0.8].

RVE determination
This section aims to determine the size for the representative volume element (RVE) in the present configuration by performing triaxial tests with different number of particles: N part ∈ {500 × i; i ∈ 1, 20 }. A reference triaxial test obtained for N max part = 20,000 will be used to compute the root mean square error (RMS error) of each simulation. For a quantity s, the RMS error s err over all the N pts points of the shearing phase is computed as follows: (4) This serie will be called S1. All the triaxial tests in this serie are performed with P c = 100 kPa. A target initial void ratio is set by slowly reducing the inter-particle friction angle during the isotropic compaction (as in [11]), all samples thus start the shearing phase with e 0 ≈ 0.60. Figure 2 shows q and V during all the simulations of S1. The deviatoric stress q is defined as the difference between the stress on the shearing axis and the lateral stress.  The volumetric strain V is defined as the trace of the true strain tensor. The RMS error based on Equation (4) was plotted for q and V on Figure 3. One can see that q err and err V are decreasing for N part ≤ 7,000 and then progressively reach a plateau. In addition, Figure 4 shows the time cost T c of all the simulations in S1, being run sequentially on a processor Intel(R) Xeon(R) CPU E5-2623 v3 @ 3.00GHz. Obviously T c increases almost linearly with N part . Considering the results given in Figure 3 and Figure 4, the best compromise between precision and time cost is met for N part = 7,500.

Critical state and rattlers' influence
Using now 7,500 particles, another serie of triaxial tests (called S2) was performed using 7 different confining pressures P c and 15 different initial void ratios e 0 , constituting a large collection of 105 triaxial tests to plot the CSL. As already explained, the initial density was set by modifying the inter-particle friction angle during the isotropic compaction. In Table 2 the values used for P c are given alongside the minimum and maximum initial void ratio obtained for each P c , namely e min 0 and e max 0 . Figure 5 shows the e 0 obtained for each ϕ comp under all P c .
The red lines on Figure 6, fitted to the critical states for P c < 40 kPa and P c ≥ 40 kPa, represent the obtained    CSL e crit (p). One can see that for P c < 40 kPa the CSL is increasing, such a counter-intuitive behaviour is also observed in [6] for ϕ > 26 • . The inter-particle friction angle indeed helps stabilising the force chains and makes some particles unnecessary to maintain the static equilibrium of the packing. Such a passive particle is called a rattler and is considered to be so if it has only 0 or 1 contact. Their number, N rattlers , is depicted in Figure 7. It can be stated that the lower P c , the higher and noisier N rattlers , and that a critical (e 0 -independent) value is also eventually obtained. Beside, N rattlers stabilises faster for loose samples but slower for high P c , except if P c is high enough for rattlers to be too rare (P c > 100 kPa). Figure 8 shows the critical mean coordination number Z crit c against the critical number of rattlers N crit rattlers normalised by the number of particles. It is noticeable that Z crit c is strongly correlated with N crit rattlers and decreases linearly for P c ≤ 40 kPa.  The great proportion of rattlers at low P c (around 20%) shows that the CSL is biased by these passive particles and suggests that they could be preferably counted as voids. To this end, a void ratio "without rattlers" e WOR was computed following [6]. Denoting V void the volume of all the voids in the sample, V part the volume of all the particles in the samples and V clp the volume of one particle, it comes: The blue lines plotted on fig. 6 corresponds to the CSL obtained by considering e WOR crit instead of e crit . Obviously, this consideration is enough to obtain a classical decreasing CSL at low P c .

Conclusion
In this paper, triaxial tests were performed on tetrahedral particles using a visco-elastic contact model with friction and no gravity. It was first determined that 7,500 particles is an appropriate sample size to constitute a RVE for the present triaxial analysis. Then, the critical states reached with triaxial tests performed at several confining pressures and initial packing densities showed that rattlers should be considered as void, specially for confining pressures lower than 40 kPa. Indeed, including rattlers into the computation of the critical packing density makes the critical state line to unexpectedly become increasing at low confining pressures. One could investigate further in the matter by performing triaxial tests under gravity, which should in principe reduce substantially the number of rattlers at low pressure. It is worth noting that the material studied here only contains identical particles, the results might differ for more classical granular materials. Indeed small and coarse particles, usually present in granular materials, should increase the number of rattlers even for high confining pressures.