Stochastic model for the micromechanics of jammed granular materials : experimental studies and numerical simulations

To describe statistical properties of complicated restructuring of the force network under isotropic compression, we measure the conditional probability distributions (CPDs) of changes of overlaps and gaps between neighboring particles by experiments and numerical simulations. We find that the CPDs obtained from experiments on a 2D granular material exhibit non-Gaussian behavior, which indicates strong correlations between the configurations of overlaps and gaps before and after deformation. We also observe the non-Gaussian CPDs by molecular dynamics simulations of frictional disks in two dimensions. In addition, the numerically calculated CPDs are well described by q-Gaussian distribution functions, where the q-indices agree well with those from experiments.


Introduction
The mechanics of granular materials has been widely investigated because of their importance in industry and science.However, mechanical properties of jammed granular materials are still not fully understood due to the disordered configurations of the constituent particles [1].At the microscopic scale, the mechanical response to quasistatic deformations manifests as changes of the contact and force network [2], where complicated rearrangements of the constituents cause irreversible restructuring.This evolution under quasi-static deformation is described by the change of the probability distribution function (PDF) of these forces, where many theoretical studies [3,4] have been devoted to determine their functional forms observed in experiments [5] and numerical simulations [6].Macroscopic quantities are then defined as statistical averages over the contact forces.
Recently, we have proposed a stochastic approach to the irreversible restructuring of the force network in bidisperse frictionless [7] and polydisperse frictional [8] particles, where the change of the PDF of both forces and interparticle "gaps" was predicted by a master equation.The master equation encompassed the mechanical response to compression (or decompression), where transition rates, or conditional probability distributions (CPDs), fully determine the statistics of the micromechanic evolution of jammed granular materials.
In this paper, we measure the CPDs for the master equation by experiments on wooden cylinders [9] and compare our experimental results with those of numerical simulations.

Microscopic responses of overlaps
When granular materials are deformed, the force network is restructured because of complicated rearrangements of the constituent particles (i.e., non-affine deformation), see Fig. 1.To detect this restructuring of the contact and force network (including opening and closing contacts), we introduce the Delaunay triangulation (DT) (Fig. 2), where not only the particles in contact, but also the nearest neighbors without contact, i.e., virtual contact, are connected by Delaunay edges.Then, we define a generalized overlap between the particles, i and j, as where R i (R j ) and D i j represent the radius of i-th ( j-th) particle and the length of the Delaunay edge connecting the particles, i and j, respectively.Note that the generalized overlap is positive if the particles are in contact, while it becomes negative if the particles are not in contact (hereinafter, virtual contact), see Fig. 1.
The affine response of generalized overlap to isotropic compression is given, at first order, by x affine i j = x i j + (D i j δφ/2φ), where the area fraction of constituent particles increases from φ to φ + δφ with the (small) increment, where the red and blue solid lines represent contact status (red for grains in contact, blue for virtual contacts).
The four kinds of transitions, (CC) contact-to-contact, (VV) virtual-to-virtual, (CV) contact-to-virtual, and (VC) virtual-to-contact, are displayed.
δφ.In granular materials, however, constituent particles are randomly arranged and force balance is violated by compression.Then, the particles exhibit a deviation from affine behaviour and complicated rearrangements such that the force balance is restated.The system relaxes to another mechanical equilibrium and the generalized overlap changes to a new value, x i j x affine i j , that is the nonaffine response of the generalized overlaps.As shown in Fig. 1, there are only four kinds of transitions from x i j to x i j : both overlaps are positive (x i j , x i j > 0) or negative (x i j , x i j < 0), if they do not change sign and thus (virtual) contacts are neither generated nor broken.We call these transitions contact-to-contact (CC) and virtual-to-virtual (VV), respectively.On the other hand, positive overlaps can change to negative (x i j > 0 and x i j < 0) and negative ones can become positive (x i j < 0 and x i j > 0), where existing contacts are broken and new contacts are generated, respectively.We call these transitions contact-to-virtual (CV) and virtual-to-contact (VC), corresponding to opening and closing contacts, respectively.
In the following, we scale the generalized overlaps by the averaged overlap before compression, x, and introduce the scaled overlaps before deformation and after relaxation as ξ ≡ x i j / x and ξ ≡ x i j / x, respectively, where we omit the subscript, i j, after the scaling.

A master equation for the PDFs of overlaps
The restructuring of the force network attributed to the transitions (CC, VV, CV, and VC) is encompassed by the evolution of the PDF of scaled overlaps, P φ (ξ), where the subscript, φ, indicates the area fraction of the system.The PDF changes to P φ+δφ (ξ ) under quasi-static isotropic compression 1 .However, it is difficult to predict the change from P φ (ξ) to P φ+δφ (ξ ) because the change from ξ to ξ is stochastic rather than deterministic. 1 The total number of Delaunay edges is conserved and the PDF is normalized as To describe this stochastic evolution of the force network (or the DT), we connect the PDF after relaxation to that before compression through the Chapman-Kolmogorov equation [10], where the transition from ξ to ξ is assumed to be a Markov process.On the right-hand-side of Eq. ( 2), the CPD of scaled overlaps has been introduced as W(ξ |ξ) 2 .Then, a master equation for the infinitesimal evolution of the PDF can be derived from Eq. ( 2) as [10]: where T (ξ |ξ) = lim δφ→0 W(ξ |ξ)/δφ has been introduced as the transition rate.The first and second terms on the right-hand-side of Eq. ( 3) represent the gain and loss of overlaps after deformation, ξ , respectively.Therefore, the transition rates or the CPDs describe the statistical properties of the micromechanic restructuring of the force network.

Experimental apparatus
In order to experimentally validate the stochastic approach presented in Sec. 2, we performed isotropic compression tests with the "1γ2 " apparatus described in [9,11].To reproduce the behaviour of a granular assembly in 2D, Schneebeli rods (i.e., 6 cm long cylindrical rollers) were used.The assembly, consisting of around 1850 rollers, is polydisperse (diameters from 8 to 20 mm).It is enclosed by a rectangular frame, see Fig. 2. Wooden rods (ash wood) were employed in order to take advantage of their quasi-elastic behaviour within the stress range applied (up to 250 kPa) at the boundaries of the sample.The grain kinematics are assessed by applying Digital Image Correlation (DIC) on a set of 80 Mpixel digital images.In particular, the technique applied, called Particle Image Tracking, is an optimisation of DIC for the case of discrete materials, for which grain motion can be erratic [9,12].
Assuming grains as 2D circles, contact deformation is modelled through the definition of a particle overlap, as in DEM, which is computed after detection of particle positions, following the definition in Eq. (1).

Numerical simulations
To model our experiments by numerical calculations, we employ molecular dynamics (MD) simulations of disks.The number of disks and the size-distribution are resembling those in our experiments, as well as the stiffness level κ, which allows us to assume that the same level of contact deformation is attained in the two cases.
The normal force between the particles in contact is modeled by a linear spring-and-dashpot, i.e., f ni j = k n y i j − 2 By definition, the CPD is normalized as η n ẏij , for y i j > 0, (having assumed a unit mass for all disks, despite the different sizes), where k n and η n are a normal stiffness and normal viscosity coefficient, respectively, and y i j is the overlap.Thus, the relative speed in the normal direction is given by its time derivative, ẏij .The tangential force is also introduced as f ti j = k t z i j − η t żij , which is switched to the sliding friction or Coulomb's friction, μ| f ni j |, when it exceeds this threshold.Here, we introduced k t = k n /2, η t = η n /4, and μ = 0.5 (consistently with the experimental value) as tangential stiffness, tangential viscosity coefficient, and friction coefficient, respectively.In addition, z i j and żij represent relative displacement and relative speed in tangential direction, respectively [13].
We first generate a dense packing of the disks in a L×L square periodic box as described in [7].Then, we apply isotropic compression to the system by rescaling every radius proportionally to the area fraction increment δφ and to the radius itself.After compression, we relax the system to mechanical equilibrium.For later analyses, we introduce a scaled strain increment as γ ≡ δφ/(φ − φ J ), with the socalled jamming density, φ J 0.8458 [14], where we have confirmed in our previous study of bidisperse frictionless particles [7] that the widths of the CPDs were proportional to γ.

Results
Fig. 3 shows experimental results of the CPDs for CC transitions, where the scaled strain increment was γ 0.074.In this figure, the CPDs are nicely symmetric and are well described by γW CC (ξ |ξ) = f c (Ξ c /γ) (the solid lines), where Ξ c ≡ ξ − m c (ξ) is the distance from the mean value, m c (ξ), and The plot refers to one single loading step.Scaling of overlaps by the strain increment, γ, was applied.The solid lines are q-Gaussian distributions with q c 1.21.
Figure 4: The CPD (a) for CC transitions with its semilogarithmic plot (b) from numerical simulations.The different symbols represent different values of the scaled strain increment, γ, as listed in the legend of (a).The solid lines are q-Gaussian distributions with q c 1.19.
is the q-Gaussian distribution function [15] with: where B(x, y) is the beta function.Then, the shape of the CPD is characterized by the q-index, q c , defined in the range between 1 < q c < 3, where we find q c 1.21 from the experiments.Therefore, the CPD for CC transitions deviates from the Gaussian distribution (corresponding to the limit q c → 1), indicating some moderate correlations between the changes of overlaps.Fig. 4 shows our numerical results of the CPDs for CC transitions, where the different symbols represent different values of the scaled strain increment, γ.In this figure, all the data well collapse if we multiply the CPDs, W CC (ξ |ξ), and distances from the mean value, Ξ c , by γ and 1/γ, respectively, which means that the CPDs are selfsimilar against the change of scaled strain increment γ, for small γ [7].The solid lines are the q-Gaussian distributions, Eq. ( 5), where the q-index is given by q c 1.19 which is very close to our experimental result.
Fig. 5 displays experimental results of the CPDs for VV transitions, where γ 0.074.The solid lines represent q-Gaussian fits, where the q-index is fitted as q v 1.31.
Fig. 6 shows our numerical calculations of the CPDs for VV transitions, where the different symbols represent Figure 5: The CPD (a) for VV transitions with its semilogarithmic plot (b) in the experimental case, as in Fig. 3.The solid lines are q-Gaussian distributions with q v 1.31.Figure 6: The CPD (a) for VV transitions with its semilogarithmic plot (b), from numerical simulations, as in Fig. 4. The solid lines are q-Gaussian distributions with q v 1.09.different values of γ.Here, we also confirm the selfsimilarity of the CPDs, where they are well collapsed by the same scaling used in Fig. 4. In this figure, the solid lines are q-Gaussian fits with the q-index, q v 1.09, which turns out to be considerably smaller than for the experiments and for frictionless data [7,8].

Summary and outlook
In this study, we have experimentally and numerically investigated the CPDs of overlaps and gaps between neighboring particles in jammed granular materials.From our experimental tests of isotropic compression, we have found that the changes of overlaps and gaps are strongly correlated such that the CPDs exhibited the non-Gaussian behavior.Assuming Tsallis' q-Gaussian statistics for the complicated restructuring, we have described the CPDs by the q-Gaussian distribution functions, where the experimentally obtained q-indices well agreed with those obtained from our MD simulations of frictional disks.
Similar q-indices are related to similar shapes of the curves in Figs.3-4; there are several possible interpretations for this similarity, as it can be related to long-range correlations (meaning that the same correlation length is found in experiments and simulations), but also to longtime correlations, memory and other alternative physical explanations.Apart from the exact interpretation, the most relevant feature to highlight is the reproducibility of some aspects between experiments and simulations.
In future, we plan to extend our study to the case of simple shear deformations, where more drastic restructuring and "recombinations" of the force network can be expected.For completeness, we also need to study the CPDs for contact-to-virtual (CV) and virtual-to-contact (VC) transitions, as in Ref. [7].

Figure 1 :
Figure 1: (Colour online) Schematic pictures of jammed packings (a) before and (b) after quasi-static deformation,where the red and blue solid lines represent contact status (red for grains in contact, blue for virtual contacts).The four kinds of transitions, (CC) contact-to-contact, (VV) virtual-to-virtual, (CV) contact-to-virtual, and (VC) virtual-to-contact, are displayed.
Colour online) Image of the initial configuration of a sample under an isotropic stress loading of 20 kPa in the "1γ2 " device.The frame is 599.4 × 444.7 mm.The solid lines represent the edges obtained with a Delaunay Triangulation: red lines are for real contacts, while virtual contacts (i.e., neighbours) are in blue.

Figure 3 :
Figure 3: The CPDs (a) for CC transitions multiplied by the scaled strain increment, γ, and plotted against scaled distances from the mean values, Ξ c /γ, with its semilogarithmic plot (b), in experiments.The plot refers to one single loading step.Scaling of overlaps by the strain increment, γ, was applied.The solid lines are q-Gaussian distributions with q c 1.21.