Understanding the effects of inter-particle contact friction on the elastic moduli of granular materials

Understanding the mechanical stiffness of closely packed, dense granular systems is of interest in many fields, such as soil mechanics, material science and physics. The main difficulty arises due to discreteness and disorder in granular materials at the microscopic scale which requires a multi-scale approach. The Discrete Element Method (DEM) is a powerful tool to inspect the influence of the microscopic contact properties of its individual constituents on the bulk behavior of granular assemblies. In this study, the isotropic deformation mode of polydisperse packings of frictionless and frictional spheres are modeled by using DEM, to investigate the effective stiffness of the granular assembly. At various volume fractions, for every sample, we determine the stress and fabric incremental response that result from the application of strain-probes. As we are interested first in the reversible, elastic response, the amplitude of the applied perturbations has to be small enough to avoid opening and closing of too many contacts, which would lead to irreversible rearrangements in the sample. Counterintuitively, with increasing inter-particle contact friction, the bulk modulus decreases systematically with the coefficient of friction for samples with the same volume fraction. We explain this by the difference in microstructure (isotropic fabric) the samples get when compressed to the same density.


INTRODUCTION
Granular materials behave differently from usual solids or fluids and show peculiar mechanical properties like dilatancy, history dependence, ratcheting and anisotropy, see [21] and references therein. The behavior of these materials is highly interesting, non-linear and involves irreversibility (plasticity), possibly already for very small strains, due to rearrangements of the elementary particles [7,17,2,11]. The concept of an initial purely elastic regime (at very small strain) for granular assemblies is an issue still under debate in the mechanical and geotechnical communities [4]. Recent works [5,22,13] show that along with the macroscopic properties (stress and volume fraction) also the structure, quantified by the fabric tensor [10,11] plays a crucial role, as it characterizes, on average, the geometric arrangement of contacts, i.e. the microstructure of the particle packing. In particular, when the material is sheared, anisotropy in the contact network develops, as related to the opening and closing of contacts, restructuring, and the creation and destruction of force-chains. The anisotropic state is at the origin of interesting observations on wave propagation in sheared granular media. However, for isotropic deformations, shear is not active so that we can focus on the interplay between contact friction and the isotropic fabric.
In order to investigate the evolution of the incremental response as a function of loading (in different directions) we performed so-called strain probing tests along an axisymmetric deformation path [12], resembling what happens in physical experiments [16]. The directionaldependent data can be collected in stress and fabric response envelopes as defined in [8] in the case of both small (elastic) and large (elasto-plastic) probes. In the case of a finite assembly of frictionless particles, in simulations, a finite elastic regime can always be detected and the elastic stiffnesses can be measured by means of an actual, very small strain perturbation [15]. In this short paper, we focus on the isotropic deformation mode only, but try to address the question on what the effects of friction are before, in ongoing research, we also address the different deformation directions axial-symmetric and plane-strain. By isolating elasticity we aim to identify the kinematics at the microscale that leads to either macroscopic elasticity or plasticity. The final goal is to improve the understanding of inter-particle contact friction on the microstructure and thus on the elastic bulk modulus in dense frictional particle systems and to guide further developments for new constitutive models based on microstructure information, as presented by [11] for frictionless samples, complementing the uni-axial compression tests by [10].

NUMERICAL SIMULATION
The Discrete Element Method (DEM) [14,11] can help to understand the deformation behavior of particle systems. At the basis of DEM are force laws that relate the interaction force to the overlap and tangential displacement of two particle contact surfaces. If all forces f i acting on particle i, from all sources, are known, the problem is reduced to the integration of Newton's equations of motion for the translational and rotational degrees of freedom, see for example [19] for more details

Contact model
For the sake of simplicity, the linear visco-elastic contact model for the normal component of force is used. The simplest normal contact force model is given by f n = kδ + γδ, where k is the spring stiffness, γ is the contact viscosity parameter, δ = (d i + d j ) /2 − (r i − r j ) .n is the overlap between two interacting species i and j with diameters d i and d j , with contact normal vector n = (r i − r j ) / |(r i − r j )| andδ is the relative velocity in the normal direction. In order to reduce dynamical effects and shorten relaxation times, an artificial viscous background dissipation force f b = −γ b v i proportional to the moving velocity v i of particle i is added, resembling the damping due to a background medium, as e.g. a fluid. The tangential force model introduced by [14] is used and thus will not be detailed here.
The standard simulation parameters are, N = 4096(= 16 3 ) particles with average radius The tangential stiffness and viscosity are set as k t /k = 0.2 and γ t /γ n = 0.2, while the coefficient of friction is varied. Note that the polydispersity of the system is quantified by the width (w = r max /r min = 3) of a uniform size distribution, where r max and r min are the radii of the biggest and smallest particles respectively. For details about other time scales present in the system, see [9,11].

Macroscopic (tensorial) quantities
Here, we define averaged tensorial macroscopic quantities -including strain-, stress-and fabric (structure) tensors -that provide information about the state of the packing and reveal interesting bulk features.
By speaking about the strain-rate tensorĖ, we refer to the external strain that we apply to the sample. The isotropic part of the infinitesimal strain tensor ε v is defined as: ε v =ε v dt = − (ε xx + ε yy + ε zz ) /3 = tr(−E)/3 = tr(−Ėdt)/3, where ε αα =ε αα dt with αα = xx, yy and zz as the diagonal components of the tensor in the Cartesian x − y − z reference system. The trace integral of 3ε v is denoted as ε v , the true or logarithmic volumetric strain, i.e., the volume change of the system, relative to the initial reference volume, V 0 . On the other hand, from DEM simulations, one can measure the 'static' stress in the system averaged over all the contacts in the volume V of the dyadic products between the contact force f c and the branch vector l c , where the contribution of the kinetic fluctuation energy has been neglected [9,10]. The isotropic component of the stress is the pressure P = tr(σ)/3. In order to characterize the geometry/structure of the static aggregate at microscopic level, we measure the fabric tensor, defined as weighted according to V P , the particle volume for particle P, which lies inside the averaging volume V , and n c is the normal unit branch-vector pointing from center of particle P to contact c [11]. Note that in this work we use k * = k/ (2 r ) to non-dimensionalize the stress, i.e. σ * = σ/k * .

Sample preparation and test procedure
In this subsection, we first describe the preparation procedure and then the details of the numerical isotropic test. The preparation procedure is an essential step in any physical/numerical experiment to obtain reproducible and reliable results, especially when friction is involved. The initial configuration is such that spherical particles are randomly generated, with low volume fraction and rather large random velocities in a periodic 3D box, such that they have sufficient space and time to exchange places and to randomize themselves. The initial configuration is obtained by, first compressing a granular gas up to a volume fraction below the jamming fraction. The system is then relaxed to allow the particles to dissipate kinetic energy and achieve a zero-pressure static configuration [9,11]. This is followed by a compression-decompression isotropic cycle to a desired maximum volume fraction ν max = 0.82 [9,11]. This procedure is applied for several packings of polydisperse spheres with different coefficients of friction while the randomization procedure was the same for all the samples. Note that the preparation is carried out with strain-control, where at every time-step the particles are first moved according to the momentary strain-rate and then allowed to relax according to the momentarily acting forces (for one time step) before the next step.

BULK MODULUS
We now study the evolution of the incremental behavior during isotropic compression as function of the micro property contact-friction. Various configurations, at different volume fraction above jamming along the unloading branch, the first one being near jamming point and last one highly dense packed, are chosen for all the packings and sufficient relaxation is applied to allow the particles to achieve a static configuration in mechanical equilibrium.
Further on we probe the samples, by applying small strain perturbations and measuring the incremental (stress and fabric) responses [15,9,10,11] to investigate the effects of inter-particle contact friction on the elastic bulk modulus. For each state we apply an identical strain-rate with amplitudeε v = 10.0 [µs −1 ] applied. After probing the configurations, the bulk modulus is obtained for the granular assembly as the ratio between the measured increment in pressure and the applied isotropic strain

Evolution of the bulk modulus
As we are interested in the elastic response, we first have to identify the elastic regime, the marginal regime and the plastic regime [11]. For the results presented below, the applied infinitesimal strain step is kept small enough to avoid strong, irreversible particle rearrangements, while plasticity develops in the sample as soon as those rearrangements happen. The response of the bulk modulus against the amplitude of the applied isotropic 3δε v strain is depicted for a chosen configuration during probing by using the so-called incremental stress response in Fig.  1 to test the rearrangements argument. It is clear that, B stays practically constant for small amplitudes (3δε v < 10 −4 ) and the regime can be considered to be elastic. By increasing the amplitudes of the small perturbation 3δε v , B starts to increase non-linearly as soon as the first change in the number of contacts happens. It is interesting to notice that for B, the elastic regime expands when the friction coefficient increases (data not shown).

Effect of inter-particle contact friction on the bulk modulus
In this subsection, the effect of inter-particle contact friction on the mechanical response of packings will be analyzed, following the strain-probes methodology presented in section 2.3. In Fig. 2(a), we plot the variation of the bulk modulus B, with volume fraction for packings with different coefficients of friction µ. For higher densities, the bulk modulus increases by decreasing the friction coefficient for samples with the same volume fraction. This is related to a higher average number of contacts for samples prepared with low friction (given the rather low tangential stiffness k t used here). For low densities, the behavior is opposite, since the higher friction at the contact leads to a lower jamming volume fraction. Thus, frictional and frictionless curves show qualitatively different slopes, associated with the activation of tangential forces. In general, one can observe three regimes, frictionless, low and high friction. However, when the data are plotted not against volume fraction, but against the isotropic fabric, see Fig.2(b), all data have about the same slope. The low friction regime behaves similar to the frictionless one at the loose states (close to the jamming volume fraction) and it follows the behavior of high friction systems when far from the jamming point.

CONCLUSIONS
In a triaxial box, the bulk modulus used to describe the incremental, elastic constitutive behavior of granular material can be measured by applying small strain perturbations to relaxed states. Those states can have experienced different deformation history, or the particles can have different properties, like the coefficient of inter-particle contact friction. The goal of this paper was to understand the effect of friction on the macroscopic bulk modulus and its relation to the micro-stucture of the packing.
Concerning the probing, the amplitude of the applied perturbations has to be small enough to avoid particle rearrangements and thus to get the elastic response, whereas large amplitudes develop plasticity in the sample due to contact and structure rearrangements between particles. A relation between the macroscopic elastic response and the micro-structure was investigated, by considering both the bulk modulus B and volumetric fabric F v , respectively. Surprisingly, frictionless packing do not follow the same trend as frictional packings. The tangential force plays a crucial role in establishing the contact network and thus for the mechanical behavior of granular materials. This is a big step towards realistic materials; since frictionless particles are an unrealistic limit-case material, the tangential force has to be taken into account, even though it carries relatively small magnitude of force, in the case of low friction coefficients.
Extension of the work to investigate the influence of inter-particle contact friction on the shear modulus is in progress, using the same small perturbation probing and allowing for many different modes of deformation. Additional future work will focus on establishing a constitutive model for frictional particles, involving the elastic moduli measured and the relation between effective moduli and the microstructure tensor (fabric).