Multi-scaleA
Master equation
Materialfor
Modelsforce distributions in polydisperse frictional particles

IV International Conference on Particle-based Methods – Fundamentals and Applications
PARTICLES 2015
E. Oñate, M. Bischoff, D.R.J. Owen, P. Wriggers & T. Zohdi (Eds)

A MASTER EQUATION FOR FORCE DISTRIBUTIONS IN
POLYDISPERSE FRICTIONAL PARTICLES
Kuniyasu Saitoh1 , Vanessa Magnanimo1 , and Stefan Luding1
1 Faculty of Engineering Technology, MESA+, University of Twente

Drienerlolaan 5, 7522 NB, Enschede, The Netherlands

Key words: Friction, Granular Materials, Force Chains, Quasi-static Deformations,
Stochastic Model, DEM
Abstract. An incremental evolution equation, i.e. a Master equation in statistical mechanics, is introduced for force distributions in polydisperse frictional particle packings.
As basic ingredients of the Master equation, the conditional probability distributions ofparticle overlaps are determined by molecular dynamics simulations. Interestingly, tailsof the distributions become much narrower in the case of frictional particles than frictionless particles, implying that correlations of overlaps are strongly reduced by microscopic
friction. Comparing different size distributions, we find that the tails are wider for thewider distribution.

1

INTRODUCTION

Quasi-static deformations of soft particles, e.g. glasses, colloids, and granular materials, have been widely investigated because of their importance in industry and science.
However, their macroscopic behaviors are still not fully understood due to disordered
configurations and complex dynamics [1]. At the microscopic scale, their mechanicalresponse is probed as a reconstruction of the force-chain network [2], where non-affine
displacements of particles cause a complicated restructuring of the network including also
opening and closing contacts. If a macroscopic quantity, e.g. stress tensor, is defined as astatistical average in the force-chain, its response to quasi-static deformations is governedby the change of probability distribution function (PDF) of forces. Therefore, the PDFs
in soft particle packings have practical importance such that many theoretical studies[3, 4] have been devoted to determine their functional forms observed in experiments [5]and simulations [6].
Recently, we have proposed a Master equation for the PDFs as a stochastic approach
towards a microscopic theory for quasi-static deformations of two-dimensional bidisperse
frictionless particles [7]. The Master equation can reproduce stochastic evolution of the
PDFs under isotropic (de)compressions, where the conditional probability distributions(CPDs) for the Master equation fully encompass the statistics of microscopic changes
11028

Kuniyasu Saitoh, Vanessa Magnanimo and Stefan Ludingof force-chain networks. In addition, any changes of macroscopic quantities constructedfrom the moments of forces can be predicted by the Master equation.
In this paper, we generalize our stochastic approach towards wider four-disperse size
distributions and include friction in the contacts. We determine the CPDs from MD
simulations and compare them with the cases of polydisperse frictionless and bidisperse
frictionless particles. First, we explain our MD simulations in Sec. 2. Then, we show ourresults in Sec. 3 and conclude in Sec. 4.
2

METHODWe use MD simulations of two-dimensional four-disperse mixtures of frictional soft
particles. The number of constituents (N1 , N2 , N3 , and N4 ) and the size distribution are
the same with those used in the experiments of wooden cylinders [8] as listed in Table
1, where the mass, m, is identical in simulations to be used for the unit of mass. The
normal force between the particles in contact is given by fijn = kn xij − ηn ẋij with a normal
stiffness, kn , and normal viscosity coefficient, ηn . Here, an overlap between the particles
i and j is introduced as
xij = Ri + Rj − dij(1)with an interparticle distance, dij , and the particles’ radii, Ri and Rj . Thus, the relative
speed in the normal direction is given by its time derivative, ẋij . The tangential force isintroduced as fijt = kt yij − ηt ẏij which is switched to the sliding friction, fijs = −µ|fijn |,when the tangential force exceeds the critical value, i.e. if |fijt | > µ|fijn |, where kt = kn /2,
ηt = ηn /4, and µ are a tangential stiffness, tangential viscosity coefficient, and friction
coefficient, respectively. Here, yij and ẏij are a relative displacement and speed in thetangential direction, respectively (see the details in Ref. [9]).
To make static packings of N (= N1 + N2 + N3 + N4 ) = 1872 particles, we randomly distribute them in a L × L square periodic box, where no particle touches others
and the friction coefficient is set to zero, µ = 0 (i.e. we use frictionless particles during the preparation of static packings). Then, we rescale every radius as Ri (t + δt) =[1 + {x̄ − xm (t)} /λ] Ri (t) (i = 1, . . . , N ), where t, δt, x̄, and xm (t) are time, an increment of time, a target value of averaged overlap, and the averaged overlap at time t,
respectively. Here, we use a long length scale λ = 102 σ̄ with the mean diameter inthe final state, σ̄, to rescale each radius gently 1 . During the rescaling, each radius increases if the averaged overlap is smaller than the target value, xm (t) < x̄, and viceversa, so that the averaged overlap will finally converge to x̄. Note that neither particle masses nor the ratios between different diameters change during the rescaling, i.e.
σi (t + δt)/σj (t + δt) = σi (t)/σj (t). We stop the rescaling when every acceleration ofparticles drops below a threshold, 10−6 kn σ̄/m, and assume that the system is static.
1

We confirmed that static packings prepared with longer length scales, λ = 103 σ̄ and 104 σ̄, give thesame results concerning their power-law behaviors, i.e. critical scaling, near jamming [10], while we cannot
obtain the same results with a shorter length scale, λ = 10σ̄.

2

1029

Kuniyasu Saitoh, Vanessa Magnanimo and Stefan Luding

Figure 1(a) is a snapshot of our simulation, where the system is static with the averaged
overlap, x̄ = 3 × 10−6 σ̄, and each color corresponds to each constituent as listed in Table
1. Figure 1(b) shows the complete Delaunay triangulation (DT) of the same packing in
Fig. 1(a), where the red solid lines represent force-chain networks, while the blue solidlines connect the nearest neighbors without contacts. In our simulations, the distancesfrom jamming are determined by the critical scaling of averaged overlap, x̄ = A(φ − φJ )[10], where the critical amplitude is found to be A ≃ 0.25σ̄ as shown in Fig. 2(a). We
also confirm other critical scalings of the static pressure divided by the normal stiffness(Fig. 2(b)) and the first peak value in radial distribution function of scaled distance (Fig.
2(c)), where we find that the power law dependence on the distance from jamming are
given by p/kn = 1.25 × (φ − φJ )1.04 and g1 = 0.31(φ − φJ )−1 , respectively 2 [10].
Then, we switch the friction coefficient to µ = 0.5 and apply
 an isotropic compression
to the prepared packings by multiplying every diameter by 1 + δφ/φ such that the area
fraction increases from φ to φ + δφ. After compression, we relax the system until every
acceleration of particles drops below the threshold again.
Table 1: The number of particles per species, Ns , diameters, σs , and masses, ms , of the four kinds ofparticles (s = 1, 2, 3, 4) resembling values used in experiments with wooden cylinders in Ref. [8]. Eachcolor in Fig. 1 are also listed in the last column.

s
1
2
3
4

3

Ns
807
434
414
217

σs /σ4
0.3
0.6
0.8
1.0

ms /m4
1.0
1.0
1.0
1.0

colors
green
red
blue
gray

RESULTSIn this section, we introduce a Master equation for the PDFs of forces as a stochastic
description of microscopic changes of force-chain networks. First, we study microscopic
responses of force-chain networks to isotropic compressions (Sec. 3.1), where we describemean values and fluctuations of overlaps in terms of applied strain increments and distances from jamming (Sec. 3.2). We then introduce a Master equation (Sec. 3.3) anddetermine transition rates for the Master equation (Sec. 3.4).
3.1

Microscopic response

At the microscopic scale in soft particle packings, the mechanical response to quasistatic deformations is probed as reconstruction of force-chain networks, where complicated2Here, p/kn is dimensionless in two-dimension and the scaled distances in radial distribution functions
are defined as r ≡ dij /(Ri + Rj ).

3

1030

Kuniyasu Saitoh, Vanessa Magnanimo and Stefan Luding
(a)

(b)

Figure 1: (Color online) (a) A static particle packing with an averaged overlap, x̄ = 3 × 10−6 σ̄, where
the total number of particles is N = 1872. (b) The complete Delaunay triangulation (DT) of the static
packing in (a): The red solid lines are equivalent to force-chains, where their widths are proportional tothe magnitudes of contact forces. The blue solid lines connect the nearest neighbors without contacts.
The gray circles represent the particles.

particle rearrangements cause the recombination of force-chains including also opening andclosing contacts. To take into account such changes in structure, we employ the complete
Delaunay triangulation (DT) of static packings as shown in Fig. 1(b), where not only theparticles in contacts, but also the nearest neighbors without contacts, i.e. virtual contacts,
are connected by Delaunay edges. We then generalize the definition of “overlaps” from
Eq. (1) toxij ≡ Ri + Rj − Dij ,(2)

where Dij is the Delaunay edge length and the overlaps between particles in virtual
contacts (Ri + Rj < Dij ) are defined as negative values. Because the DT is unique foreach packing, contacts and virtual contacts are uniquely determined.
If we apply an isotropic affine compression to the system, every generalized overlap,
Eq. (2), (not only contacts, but also virtual contacts) changes toxaffine= xij +
ij

Dij
δφ ,
2φ

(3)

where we neglected the higher order terms proportional to xij δφ and δφ2 . However,the particles are randomly arranged and the force balance is broken for each particle by
the affine deformation. Therefore, the system relaxes to a new static state, where nonaffine displacements of particles cause complex changes of contacts including opening andclosing contacts. After the relaxation, overlaps change to new values, x′ij �= xaffine, that isij
non-affine responses of overlaps.
41031

Kuniyasu Saitoh, Vanessa Magnanimo and Stefan Luding(a) 100(b) 100(c) 10

-1

10

-1

10

-2

10

-3

10

-4

10

-5

10

10

-2

10

10

-4

10

10
10

-5

10
10

5
4

10

-3

-6

6

3
2
1

10
-6

10

-5

10

-4

10

-3

10

-2

10

-1

10

0

10

-6

10

-6

10

-5

10

-4

10

-3

10

-2

10

-1

10

0

10

0

10

-6

10

-5

10

-4

10

-3

10

-2

10

-1

10

0

10

Figure 2: (Color online) Double logarithmic plots of the (a) averaged overlap, x̄, (b) static pressure
scaled by the normal stiffness, p/kn , and (c) first peak of radial distribution function for scaled distance,
g1 , where the dotted lines are power law fittings, (a) x̄ = 0.25σ̄(φ − φJ ), (b) p/kn = 1.25 × (φ − φJ )1.04 ,and (c) g1 = 0.31(φ − φJ )−1 , respectively.
(a)

(b)

(VC)

(c)

(CC)

(d)
6

(VC)

(CC)

(VV)

(CV)

3
0

(VC)

(CC)

(VV)

(CV)

-3
-6

(VV)

(CV)

-6

-3

0

3

6

Figure 3: (Color online) (a) and (b): Sketches of static packings (a) before and (b) after (de)compression,
where the red and blue solid lines represent contacts and virtual contacts, respectively. The four kinds oftransitions, (CC) contact-to-contact, (VV) virtual-to-virtual, (CV) contact-to-virtual, and (VC) virtualto-contact, are displayed. (c) A scatter plot of scaled overlaps, where the red and blue dots are ξ ′ andξ affine plotted against ξ, respectively. (d) A schematic picture of affine and non-affine responses of scaled
overlaps, where the blue and red solid lines represent ξ affine (in average) and the linear functions, ms (ξ),
respectively. The slope in (CC), ac , and all the dimensionless lengths, bs and vs , are proportional to γ,
while the slope in (VV) is negligible, av ≃ 0.

As shown in Figs. 3(a) and (b), there are only four kinds of transitions from xij to x′ij : A
positive overlap, xij > 0, remains as positive, x′ij > 0, or a negative overlap, xij < 0, stays
in negative, x′ij < 0, where they do not change their signs and thus contacts are neither
generated nor broken. We call these changes “contact-to-contact (CC)” and “virtual-tovirtual (VV)”, respectively. The other cases are that a positive overlap changes to anegative one and a negative overlap becomes positive, where existing contacts are broken
and new contacts are generated, respectively. We call these changes opening and closing
contacts, or in analogy to the above, “contact-to-virtual (CV)” and “virtual-to-contact(VC)”, respectively.
In the following, we scale the generalized overlaps by the averaged overlap before compression such that scaled overlaps before compression, after affine deformation, and after
relaxation are introduced as ξ ≡ xij /x̄, ξ affine ≡ xaffine/x̄, and ξ ′ ≡ x′ij /x̄, respectively
ij
5

1032

Kuniyasu Saitoh, Vanessa Magnanimo and Stefan Luding(we omit the subscript, ij, after the scaling). From Eq. (3) and the critical scaling,x̄ = A(φ − φJ ), the scaled overlap after affine deformation is found to beξ affine = ξ + Ba γ

(4)which is a linear function of ξ with an offset, Ba γ ≡ (Dij /2Aφ) γ, proportional to thescaled strain increment, γ ≡ δφ/ (φ − φJ ). However, the scaled overlap after non-affine deformation, ξ ′ , fluctuates around mean due to complicated particle rearrangements during
the relaxation.
3.2

Mean and fluctuation

To describe scaled overlaps after non-affine deformation, ξ ′ , we measure their mean
and fluctuations through scatter plots. Figure 3(c) displays the scatter plot, where thefour different transitions under compression are mapped onto four regions: (CC) ξ, ξ ′ > 0,(VV) ξ, ξ ′ < 0, (CV) ξ > 0, ξ ′ < 0, and (VC) ξ < 0, ξ ′ > 0, respectively. In this figure,overlaps after affine deformation are described by the deterministic equation (4), while
those after non-affine deformation distribute around mean values with finite fluctuations.
The differences between affine and non-affine responses are always present, but not visible
if the applied strain is small or the system is far from jamming, i.e. if γ ≪ 1, while ξ ′deviates more from ξ affine and data points are more dispersed if γ ≫ 1.
In the same way as affine response, Eq. (4), we describe the mean values of ξ ′ in (CC)and (VV) regions by linear functions of ξ,
ms (ξ) = (as + 1)ξ + bs ,(5)

where the subscripts, s = c and v, represent the mean values in (CC) and (VV), respectively. We also introduce standard deviations of ξ ′ from their means as vs (which
are almost independent of ξ). Then, the systematic deviation from the affine response isquantified by the coefficients, as , bs , and vs , as summarized in Fig. 3(d). Note that theaffine response, Eq. (4), is obtained if as = vs = 0 and bs = Ba γ. Except for av ≃ 0, all
the coefficients linearly increase with the scaled strain increment (Fig. 4), where all datawith a wide variety of δφ and φ − φJ collapse onto linear scalings,
as = As γ , bs = Bs γ , vs = Vs γ ,(6)between the fitting range, 10−6 ≤ γ ≤ 5 × 10−3 , with the scaling amplitudes, As , Bs ,and Vs , listed in Table 2. We observe that av ≃ 0 and Bv ≈ Ba (≃ 1.3 in average) suchthat virtual contacts almost behave affine in average except for their huge fluctuations(Vv ≫ Vc ). In contrast, Bc is always smaller than Ba such that mc (ξ) intersects ξ affine
at ξ ∗ = (Ba − Bc )/Ac ≃ 1.3, which is independent of γ. This leads to small responses,
ξ ′ < ξaffine , of small overlaps, ξ < ξ ∗ , and vice versa, implying preferred tangential andhindered normal displacements as a sign of non-affine deformations [11].
61033

Kuniyasu Saitoh, Vanessa Magnanimo and Stefan Luding
102
101

102

102
101

(a)

101

(b)

(c)

100

100

100

10-1

10-1

10-1

10-2

10-2

10-2

10-3

10-3

10-3

10-4

10-4

10-4

10-5

10-5

10-6
10-610-510-410-310-2 10-1 100 101 102

10-6
10-610-510-410-310-2 10-1 100 101 102

102

102

4x10-44x10-54x10-610-5

10-6
10-610-510-410-310-2 10-1 100 101 102
1

101

(d)
0.5
0
-0.5

100

100

10-1

10-1

10-2

10-2

10-3

10-3

-4

10-4

10-5

10-5

10

-1
10-710-610-510-410-310-210-1100 101 102

101

(e)

10

-6

10-610-510-410-310-2 10-1 100 101 102

(f)

10-6
10-610-510-410-310-2 10-1 100 101 102

Figure 4: (Color online) Double logarithmic plots of the coefficients for mean values and fluctuations ofscaled overlaps as functions of the scaled strain increment, (a) ac , (b) bc , (c) vc , (e) bv , and (f) vv , and asemi-logarithmic plot of (d) av , where we apply strain increments δφ = 4 × 10−4 , 4 × 10−5 , and 4 × 10−6
(as indicated by the different symbols in the legend of (a)) to the static packings with φ−φJ = 1.2×10−1 ,4 × 10−2 , 1.2 × 10−2 , 4 × 10−3 , 1.2 × 10−3 , and 4 × 10−4 . The solid lines represent the linear scalings,
Eq. (6), between the fitting range, 10−6 ≤ γ ≤ 5 × 10−3 .

In contrast to (CC) and (VV), the data of ξ ′ in (VC) and (CV) are concentrated in
narrow regions (the inside of the dashed lines in Fig. 3(d)), whereas ξ affine linearly increases
with ξ in (VC) and there is no data of ξ affine in (CV), i.e. affine responses do not generate
opening contacts.
Table 2: Scaling amplitudes in Eq. (6), As , Bs , and Vs , for polydisperse frictional particles. The q-indicesfor the CPDs in (CC) and (VV), where qsBL , qsPL , and qsPN are the q-indices for bidisperse frictionlessparticles [7], polydisperse frictionless particles, and polydisperse frictional particles, respectively.

s
c
v

3.3

As
Bs
0.39 0.81
0.00 1.32

Vs
qsBL
0.30 1.13
1.22 1.39

qsPL
1.72
1.96

qsPN
1.19
1.09

Master equation

The reconstruction of force-chain networks attributed to the changes, (CC), (VV),(CV), and (VC), is well captured by the PDFs of scaled overlaps. Here, we introduce
the PDF as Pφ (ξ) with the subscript, φ, representing the area fraction in the system.
Because the total number of Delaunay edges is conserved during deformations, the PDFs
71034

Kuniyasu Saitoh, Vanessa Magnanimo and Stefan Luding

0.5
0.4

PDF

0.02

0.3
0.01

0.2

0

0.1
0

-6

-6 -5 -4 -3 -2 -1 0

-4

-2

0

2

4

6

scaled overlap

Figure 5: (Color online) The PDFs of scaled overlaps before compression, Pφ (ξ) (the squares), after
affine deformation, Pφ+δφ (ξ affine ) (the triangles), and after non-affine deformation, Pφ+δφ (ξ ′ ) (the circles),where the inset is a magnification of the PDFs of negative scaled overlaps.

∞are normalized as −∞ Pφ (ξ)dξ = 1.
In our previous study of bidisperse frictionless particles [7], we found that the affine
deformation, Eq. (4), just shifted the PDF to the positive direction, while the non-affine
deformation broadened the PDF in positive overlaps and generated a discontinuous “gap”
around zero 3 . We also find that the PDF of negative overlaps after non-affine deformation is comparable with that after affine deformation. In our simulations of polydisperse
frictional particles (µ = 0.5), we observe similar results to our previous study. Figure 5
displays the PDFs for polydisperse frictional particles, where the PDF before compression, Pφ (ξ), has a discontinuous gap around zero. The affine deformation pushes the
PDF towards the positive direction as Pφ+δφ (ξ affine ), where the discontinuity is smoothed
out because of the polydispersity in the system. After non-affine deformation, the PDF
widens in positive overlaps as Pφ+δφ (ξ ′ ), while there is no significant difference between
Pφ+δφ (ξ affine ) and Pφ+δφ (ξ ′ ) in negative overlaps (the inset in Fig. 5).
To describe such the non-affine evolution of the PDFs, we connect the PDF after nonaffine deformation to that before compression through the Chapman-Kolmogorov equation
[13],

∞

Pφ+δφ (ξ ′ ) =−∞W (ξ ′ |ξ)Pφ (ξ)dξ ,(7)

assuming that transitions between overlaps (from ξ to ξ ′ ) can be regarded as Markov processes. On the right-hand-side of the Chapman-Kolmogorov equation (7), the conditional
probability distribution (CPD) of scaled overlaps, ξ ′ , which were ξ before compression,∞is introduced as W (ξ ′ |ξ). By definition, the CPD is normalized as −∞ W (ξ ′ |ξ)dξ ′ = 1.
Then, the Master equation for the PDFs is derived from Eq. (7) as [13]3Such a discontinuity is specific to static packings, where a corresponding gap has been observed in aradial distribution function of glass with zero-temperature [12].

8

1035

Kuniyasu Saitoh, Vanessa Magnanimo and Stefan Luding∂
Pφ (ξ ′ ) =
∂φ

 ∞

−∞

{T (ξ ′ |ξ)Pφ (ξ) − T (ξ|ξ ′ )Pφ (ξ ′ )} dξ(8)with the transition rate defined as T (ξ ′ |ξ) = limδφ→0 W (ξ ′ |ξ)/δφ. The first and second
terms on the right-hand-side of the Master equation (8) represent the gain and loss of new
overlaps, ξ ′ , respectively. Therefore, the transition rates or the CPDs fully determine thestatistics of microscopic responses of force-chain networks.
3.4

Conditional probability distributions

We determine the CPDs of scaled overlaps as the distributions of ξ ′ around their meanvalues, ms (ξ). For example, the CPD for affine deformation is given by a delta function,
Waffine (ξ ′ |ξ) = δ(ξ ′ −ξ affine ), which just shifts the PDF by Ba γ, i.e. Pφ+δφ (ξ) = Pφ (ξ −Ba γ).
However, non-affine deformations generate fluctuations of scaled overlaps around theirmean values so that the CPDs must have finite widths. In the following, we determine the
CPDs for polydisperse frictionless particles and those for polydisperse frictional particles
to examine the effect of particle friction (characterized by the friction coefficient, µ) on thestatistics of microscopic responses of force-chain networks. We also compare our resultswith our previous work of bidisperse frictionless particles [7] to study the influence of size
distributions.
Figures 6(a) and (c) show the CPDs for polydisperse frictionless particles (µ = 0) in(CC) and (VV), respectively. As can be seen, all data are symmetric around the meanvalues, ms (ξ), and are well collapsed if we multiply the CPDs and the distances fromthe mean values defined as Ξs ≡ ξ ′ − ms (ξ) by γ and 1/γ, respectively. In these figures,the solid lines are given by γWCC (ξ ′ |ξ) = fc (Ξc /γ) and γWV V (ξ ′ |ξ) = fv (Ξv /γ) with theq-Gaussian distribution [14],
1

 1−q
s
x2
1
fs (x) =1+,2c(qs )n(qs )Vs(9)
where the functions are defined as n(t) = (t−3)/(1−t) and c(t) = Vs n(t)B (1/2, n(t)/2)
with the beta function, B(x, y). In Table 2, we list the q-indices which characterize shapesof the CPDs. Note that the q-index must be in the range between 1 < qs < 3, where
the normal (Gaussian) distribution corresponds to the limit, q → 1. Figures 6(b) and(d) are semi-logarithmic plots of Fig. 6(a) and (c), respectively, where the dotted lines
are the CPDs obtained from our previous study on bidisperse frictionless particles [7].
From these results, we observe that the basic properties of CPDs, i.e. their symmetry andself-similarity for different γ, are not affected by the size distribution. However, we findthat the CPDs for polydisperse particles have much wider tails than those for bidisperse
particles, i.e. qsBL ≪ qsPL in Table 2, implying that the spatial correlation of scaled overlaps
increases with the increase of polydispersity [15].
Figures 7(a) and (c) display the CPDs for polydisperse frictional particles (µ = 0.5
with the same polydispersity as in Fig. 6) for (CC) and (VV) cases, respectively. Note
91036

Kuniyasu Saitoh, Vanessa Magnanimo and Stefan Luding
101

1.5

(a)

(b)
100

1

10-1
0.5

0

10-2

-2

-1

0

1

10-3

2

-2

-1

0

1

2

-10

0

10

20

100

0.15

(c)

(d)
-1

10

0.1

10-2
0.05

0
-20

10-3

-10

0

10

20

10-4
-20

Figure 6: (Color online) The CPDs for polydisperse frictionless particles multiplied by γ in (a) (CC)and (c) (VV) plotted against scaled distances from mean values, Ξs /γ (s = c, v). (b) and (d) are semilogarithmic plots of (a) and (c), respectively. The solid lines are the q-Gaussian distributions with theq-indices, (b) qcPL = 1.19 and (d) qvPL = 1.09, respectively. The dotted lines are the CPDs obtained fromour previous study of two-dimensional bidisperse frictionless particles [7], where the q-indices are givenby (b) qcBL = 1.13 and (d) qvBL = 1.39, respectively.

that all data are symmetric around their mean values and well collapse after the samescaling as in Fig. 6. Figures 7(b) and (d) are semi-logarithmic plots of Fig. 7(a) and(c), respectively, where the dotted lines are the CPDs for bidisperse frictionless particles[7]. Comparing the results in Figs. 6 and 7, we confirm that the basic properties, i.e.
symmetry and self-similarity, of the CPDs do not change by particle friction, while thetails of CPDs for frictional particles are much narrower than those for frictionless particles,
i.e. qsPN ≪ qsPL in Table 2. Therefore, the microscopic friction at the contact drasticallydecrease the spatial correlations of scaled overlaps.
4

SUMMARY

We have introduced a Master equation for the PDFs of forces in two-dimensional
polydisperse frictionless and frictional particle packings. We find that the microscopic
response to isotropic compression increasingly deviates from the affine prediction with
the increase of the scaled strain increment, γ, in agreement with our previous study ofbidisperse frictionless particles [7]. The scaling amplitudes for mean values of scaled
overlaps are almost the same as those for bidisperse frictionless particles [7], implying
10

1037

Kuniyasu Saitoh, Vanessa Magnanimo and Stefan Luding

1.5

101

4x10-44x10-54x10-6

(a)

(b)
100

1

10-1
0.5

0

10-2

-2

-1

0

1

10-3

2

-2

-1

0

1

2

-10

0

10

20

10-1

0.08

(c)

(d)
-2

0.06

10

0.04

10-3

0.02

10-4

0
-20

-10

0

10

20

10-5
-20

Figure 7: (Color online) The CPDs for polydisperse frictional particles multiplied by γ in (a) (CC)and (c) (VV) plotted against scaled distances from mean values, Ξs /γ (s = c, v). (b) and (d) are semilogarithmic plots of (a) and (c), respectively. The solid lines are the q-Gaussian distributions with theq-indices, qcPN = 1.19 and qvPN = 1.09, respectively. The dotted lines are the CPDs obtained from ourprevious study of two-dimensional bidisperse frictionless particles [7], where the q-indices are given by(b) qcBL = 1.13 and (d) qvBL = 1.39, respectively.

that the degree of non-affinity depends on neither particle friction nor size distribution
in average. In addition, symmetry and self-similarity of the CPDs, which govern thestatistics of microscopic changes of the force-chain network, are not affected by friction
and size distribution. In contrast, the shapes and widths of the CPDs are different fordifferent systems: By using the q-indices to describe their shapes in all cases, we observethat the q-index for polydisperse particles is larger than bidisperse particles, which meansthat the spatial correlation of scaled overlaps increases with polydispersity (the q-index
not equal to one means a non-Gaussian distribution). On the other hand, the q-index
becomes nearly equal to one if microscopic friction is considered, implying a dramatic
decrease of spatial correlation.
In conclusion, the basic properties of the Master equation, i.e. linear scalings of meanvalues and fluctuations of scaled overlaps, and the symmetry and self-similarity of theCPDs, are insensitive to particle friction and size distributions, in the cases tested in thispaper, while the shape of CPDs depends on both friction coefficient and polydispersity,
implying different spatial correlations of scaled overlaps for different material properties.

11

1038

Kuniyasu Saitoh, Vanessa Magnanimo and Stefan Luding

Acknowledgment
This work was financially supported by the NWO-STW VICI grant 10828. A part ofnumerical computation was carried out at the Yukawa Institute Computer Facility, Kyoto,
Japan.
REFERENCES[1] J. Lemaitre and J.-L. Chaboche, Mechanics of Solid Materials (1990). Cambridge
University Press, Cambridge, UK.
[2] T. S. Majmudar, M. Sperl, S. Luding, and R. P. Behringer, Phys. Rev. Lett. (2007)
98, 058001.
[3] S. Henkes and B. Chakraborty, Phys. Rev. E (2009) 79, 061301.
[4] J. H. Snoeijer, T. J. H. Vlugt, M. van Hecke, and W. van Saarloos, Phys. Rev. Lett.
(2004) 92, 054302.
[5] E. I. Corwin, H. M. Jaeger, and S. R. Nagel, Nature (2005) 435, 1075.
[6] L. E. Silbert, G. S. Grest, and J. W. Landry, Phys. Rev. E (2002) 66, 061303.
[7] K. Saitoh, V. Magnanimo, and S. Luding, Soft Matter (2015) 11, 1253.
[8] G. Combe, V. Richefeu, and G. Viggiani, Geomechanics from Micro to Macro, Soga
et al. (Eds). Taylor & Francis Group, London, ISBN 978-1-138-02707-7 (2015).
[9] S. Luding, J. Phys.: Condens. Matter (2005) 17, S2623.
[10] M. van Hecke, J. Phys.: Condens. Matter (2010) 22, 033101.
[11] W. G. Ellenbroek, M. van Hecke, and W. van Saarloos, Phys. Rev. E (2009) 80,
061307.
[12] L. Berthier, H. Jacquin, and F. Zamponi, Phys. Rev. E (2011) 84, 051103.
[13] N. G. van Kampen, Stochastic Processes in Physics and Chemistry, 3rd edition(2007). Elsevier B. V. Amsterdam, The Netherlands.
[14] C. Tsallis, J. Stat. Phys. (1988) 52, 479.
[15] V. Ogarko and S. Luding, Soft Matter (2013) 9, 9530.

121039