Quasistatic response of loose cohesive granular materials

DEM-simulated model cohesive assemblies of spherical grains of diameter d, with contact tensile strength F0, once prepared in loose states, are quasistatically subjected to growing isotropic pressure P, and then to triaxial compression, maintaining lateral stresses 2 = 3 = P while increasing axial stress 1 = P + q and strain ✏1. Reduced pressure P⇤ = dP/F0 varies from 0.1 (cohesion dominated case, for which systems typically equilibrate with solid fraction ' 0.35), to large values for which the cohesionless behavior is retrieved. In triaxial compression, while the moderate strain response (✏1 ⇠ 0.1) is influenced by initial coordination numbers and mesoscale heterogeneities, the approach to the critical state, as both q (deviator) and steadily increase, gets slower for smaller P⇤. Critical ratio q/P strongly increases for decreasing P⇤, as roughly predicted in an “e↵ective stress” scheme. Anomalously small elastic moduli are observed in the gel-like structures. While extensive geometric rearrangements take place, no shear banding is observed. Loose cohesive granular assemblies are thus capable of large quasistatic stable plastic strains and ductile rupture.


Introduction
Numerical grain-level investigations, by discrete element methods (DEM) [1], of the properties of assemblies of hard objects interacting by frictional contacts [2] are now quite widespread, and their quasistatic mechanical behavior, as probed, e.g., by the standard triaxial test used in geotechnique labs with sands [3], is often addressed. Cohesive granular assemblies, such as powders and colloids, are less frequently studied, but exhibit a wider variety of static states. While the solid fraction of cohesionless assemblies of spherical identical beads ranges from the random close packing value ' 0.64 down to some frictiondependent minimum, in the 0.55-0.6 range, much looser static cohesive assemblies may be observed [4], in which applied stresses are carried by tenuous, ramified contact networks reminiscent of colloidal gels [5,6].
The present contribution reports on a DEM investigation of the behavior of a simple model material, the properties of which are briefly recalled in Sec. 2, prepared in P ⇤ -dependent states by isotropic compression (Sec 3), and subjected to triaxial compression (Sec 4). A quick final discussion (Sec 5) suggests perspectives.

Particle interactions
We use the same model system as in Refs. [4,7]: spherical beads of diameter d and mass m, with Hertz contact elasticity, friction coe cient µ = 0.3, and adhesive forces chosen according to a simplified model of capillary attraction ⇤ e-mail: walid.lammali@univ-ei↵el.fr ⇤⇤ e-mail: jean-noel.roux@univ-ei↵el.fr ⇤⇤⇤ e-mail: anh-minh.tang@enpc.fr through small liquid bridges, which form in wet granular assemblies at low saturation. This results in an attractive contact force F 0 = ⇡ d ( denoting the interfacial tension of the wetting liquid). One specific feature of this force is its hysteresis: liquid bridges formed at intergranular contacts survive between receding pairs, as long as the separation distance D does not exceed a rupture threshold D R (here, D R = d/10, unless specified otherwise, corresponding to a meniscus volume 10 3 d 3 ). One should distinguish the contact coordination contribution, z c , and the distant interactions one, z d , to the total coordination number z = z c + z d . The attractive force decrease as D grows from zero to D R is modeled with a simple law [4,7] ("Maugis approximation"). The Coulomb inequality involves the elastic repulsive normal force F e N only. Thus, a contact with vanishing normal force (F e N F 0 = 0) may transmit a tangential force as large as µF 0 . Remarkably, this simple numerical model agrees quantitatively with experimental observations in simple shear flow [8,9] and in isotropic compression [4]. We use the same dimensionless control parameter as in Refs. [4,8,9], reduced pressure P ⇤ = d 2 P/F 0 , comparing the adhesive force to pressure P. Cohesive e↵ects are strong for small P ⇤ , stabilizing tenuous, open structures [4]. Their influence gradually vanish as P ⇤ increases, until the properties of cohesionless systems are retrieved for P ⇤ 1 [4,10].

Sample preparation
Our DEM procedure involve several 8000 particle samples for each set of parameters, all enclosed in boxes with fully periodic boundary conditions. To compress samples to large enough 'axial' strains ✏ 1 (counted positively) in triaxial tests, tall rectangular parallelipipedic samples are prepared, with aspect ratio 2. Loose systems need to form aggregated structures to support stresses [4,6]. Our aggregation procedure is ballistic: grains, randomly distributed within the simulation cell, at solid fraction = 0.3, are attributed Maxwell-distributed random velocities, and left to collide and interact until a unique, equilibrated cluster of aggregated grains spans the system. The initial mean quadratic velocity determines the final connectivity. F 1 = max(d 2 P, F 0 ) defining a relevant force scale, our equilibrium condition requests forces and torques to balance within respective tolerances 10 4 F 1 and 10 4 F 1 d.

Isotropic compression
Once the initial aggregated state is achieved, the sample is subjected to stepwise increasing isotropic pressure values, and equilibrated at di↵erent P ⇤ levels The resulting com- pression curve (Fig. 1) is found independent of the details of the procedure, but remains sensitive to the initial solid fraction 0 , and to the initial connectivity (see also [4]). Fig. 2 shows how better coordinated systems (right im- age), obtained with a stronger agitation in the initial aggregation stage (i.e., larger initial kinetic energy per grain, compared to F 0 D r [4,6]), comprise larger holes and dense regions. (The contact network only seems disconnected because a slab of thickness 3d is shown). The Delaunay tesselation of the set of grain centres enables the partition of the void space into disjoint pores. The decrease of pore size in the compression is shown in Fig 3, which quantifies the enduring e↵ects of the initial aggregation-induced geometric di↵erence between the states depicted in Fig. 2. Such tenuous force-carrying structures have anomalous elastic moduli, expressing their response to small stress increments without network rearrangement or irreversible sliding in contacts. This is shown in Fig. 4, in which bulk (B) and shear (G) moduli, computed with the matrix method of Refs. [10,11], are normalised by reference values B e , G e . Those are Voigt (homogeneous strain) estimates of the moduli, and their values, for Hertz-Mindlin and a material with Young modulus E and Poisson ratio ⌫, read In Eq. 1, P e denotes the elastic normal force contribution to the average stress, given by [7] P e = P + z F 0 Ratios B/B e and G/G e , anomalously small in loose states, grow with P ⇤ , and finally approach values typical of lowcoordination cohesionless grain packs [11].

Triaxial compression 4.1 Macroscopic behaviour
Trixial compression tests are applied to several states at di↵erent P ⇤ values, from 10 1 to large ones. Variations of deviator q (normalised by initial pressure P 0 ) and solid fraction , versus axial strain ✏ 1 are shown in Fig. 5. The test is carried out on controlling strain rate✏ 1 , such that the inertial number I =✏ 1 p m/dP be equal to 10 4 . It is pursued until all measured values approach the plateau corresponding to the critical state, at strains ✏ 1 which increase as P ⇤ decreases, from ⇠ 0.35 in cohesionless systems up to about 0.8 at P ⇤ = 0.1. The e↵ect of the initial coordination number (and density heterogeneities, see Figs. 2 and 3) only vanishes after this large strain interval. Interestingly, distant attractive forces have quite a notable influence on the critical shear resistance [12], but hardly a↵ect the critical density. We did not record any strain localization, in contrast with the observation of shear bands at P ⇤ = 0.1 in enduring quasistatic shear flow [7].

Critical state and yield criterion
The ratio of deviator stress q to average stress p = ( 1 + 2 3 )/3 = P 0 + q/3 would vary linearly with 1/p in a simple Mohr-Coulomb model of plastic yielding at the critical state, with cohesion c and internal friction angle ': growth of q/p with 1 , which is overestimated by the Mohr-Coulomb form identified at large p ⇤ . Yet a rough estimate of q p is obtained [7,8] on assuming a P ⇤ -independent ratio of q to the "e↵ective" average stress due to repulsive elastic forces, introduced as P e in Eq. 2: Fig. 6 shows that this prediction, although not very accurate for p ⇤ < 1, captures the correct trend.

Force networks
Coordination numbers z c and z d are shown in Fig. 7. Both approach a critical state plateau. The plateau value of z c exhibits little dependence on P ⇤ , despite the di↵erence in densities, with however a somewhat lower value in the absence of distant attraction. z d , on the other hand, increases by a larger amount in denser systems under larger P ⇤ (initial pressure and constant lateral stress). But the mechanical influence of those distant forces then decreases as 1/P ⇤ , and the contribution of distant forces to the deviatoric stress is maximal (about 10%) among simulated tests for P ⇤ = 0.1 at large strain, down to a few percents for P ⇤ 1. The major contribution (at least 75%) is given by repulsive elastic forces, with tangential forces supplying typically 20%. The influence of distant forces on the critical value of the deviator stress, as apparent in Fig. 5, is mostly indirect, probably due to the stabilization of contact networks: without distant interactions, z c decreases to the smaller value of cohesionless systems. Although we do not detail this approach here, it may be shown that the deviatoric stress is simply related to anisotropy parameters for fabric tensors and force intensities, as often pointed out for a variety of granular systems [8,13]. In this respect, the specific feature of the loose systems studied here for smaller P ⇤ values is the larger contribution of the anisotropy of forces compared to fabric terms: the contacts oriented near the major principal stress direction tend to carry larger forces, but they are not very much more numerous than those oriented differently. This might be attributed to the structures shown in Fig. 2 likely inducing correlated displacements determined by the mesoscale heterogeneities.
Elastic moduli [14,15] enable a probe of contact network anisotropy. In Fig. 8 the ratio of longitudinal moduli (those determining the velocity of longitudinal waves) in the axial and lateral directions is plotted versus ✏ 1 . Relations of elastic moduli to stress and fabric anisotropies could be exploited in experiments on such systems.

Discussion
Compared to cohesionless ones, cohesive granular assemblies exhibit a wider variety of structures and properties, and the present communication contributed to the investigation of the little known open, loose systems, with lower reduced pressure P ⇤ than in most published studies [12]. As P ⇤ increases, irreversible compaction takes place. In triaxial compression, the route to the critical state is longer, involving important density increases. The deviator to mean stress ratio achieves large values, as roughly predicted by an "e↵ective stress" approach. We showed the enduring e↵ect (i.e., up to large strains) of the initial geometry (connectivity and pore size distribution, as determined by initial aggregation process) for the same density, in isotropic or triaxial compression.
The preliminary results presented here should be further completed and analysed in many ways. The slow plastic strengthening and density increase in triaxial compression might be compared to the observed plastic response of gel networks, although those are usually probed at constant volume rather than constant normal (lateral) stress and often modeled di↵erently (e.g., with bonds and angular elasticity [16]). Elastic properties of tenuous networks, and their rupture properties call for detailed investigations. Rolling resistance in contacts [10], which strongly a↵ects loose cohesive structures, and plays an important in colloidal aggregates [17], should be introduced.
Most importantly, the mesoscale density heterogeneities of loose aggregated structures should be related to their mechanical properties.