Symposium on Large Fluctuations and Collective Phenomena in Materials III
Materials Science & Technology 2016

A MASTER EQUATION FOR FORCE DISTRIBUTIONS IN
DENSE GRANULAR MATERIALS
Kuniyasu Saitoh1 and Fumiko Ogushi1
1

WPI Advanced Institute for Materials Research, Tohoku University
2-1-1 Katahira, Aoba-ku, Sendai 980-8577, Japan

Key words: granular materials, force-chain network, molecular dynamics simulation
Abstract. Employing molecular dynamics simulations of two-dimensional soft particle packings, we investigate microscopic mechanical responses of granular materials to
simple shear deformation. Though microscopic aﬃne responses are deterministic, the restructuring of force-chain networks is well described by a stochastic method for the force
distributions, i.e. a master equation. We find that microscopic responses of force-chainsare anisotropic in average, while fluctuations of forces are equiprobable in all directions,i.e. the probability of force fluctuations is isotropic under simple shear deformation.

1

Introduction

Quasi-static deformations of granular materials have been widely investigated becauseof their importance in industry and science. However, the macroscopic behavior of dense
granular materials is still not fully understood due to disordered configurations and complex dynamics of the constituent particles [1]. At microscopic scale, mechanical responsesof granular materials are probed as the change of force-chain networks [2, 3], where complicated non-aﬃne displacements of the particles cause the restructuring of force-chains
including opening and closing contacts. If macroscopic quantities, e.g. shear stress, are defined as statistical averages in the force-chains, their mechanical responses to quasi-static
deformations are governed by the change of probability distribution function (PDF) offorces. Therefore, the PDFs in dense granular materials have practical importance suchthat many theoretical studies [4, 5] have been devoted to determine their functional forms
observed in experiments [6] and numerical simulations [7].
Recently, we have proposed a master equation for the PDFs as a stochastic descriptionof the restructuring of force-chain networks [8]. The master equation can reproduce
the evolution of the PDFs under isotropic (de)compressions, where transition rates orconditional probability distributions (CPDs) in the master equation fully encompass thestatistics of microscopic mechanics of dense granular materials.
In this paper, we apply our master equation to the stochastic evolution of force-chain
networks under simple shear deformations. First, we explain our numerical method in
1

Kuniyasu Saitoh and Fumiko Ogushi

Sec. 2. Then, we introduce our theoretical framework in Sec. 3 and show our numericalresults in Sec. 4. In Sec. 5, we discuss and conclude our results.
2

MethodWe use molecular dynamics (MD) simulations of two-dimensional 50 : 50 binary mixtures of frictionless soft particles with the same mass, m, and two kinds of radii, R and˙R/ρ (ρ = 1.4). The normal force between the particles in contact is given by f = Kξ −η ξ,˙where K, η, ξ, and ξ are a spring constant, viscosity coeﬃcient, particle overlap, and relative speed in the normal direction, respectively [9]. A global damping force, f damp = −ηv,
proportional to the particle’s velocity, v, is also introduced to enhance the relaxation,
where the damping coeﬃcient, η, is the same with the viscosity coeﬃcient between theparticles in contact. Then, we randomly distribute the N = 8192 particles in a L × L
square periodic box to make static packings by the method described in Ref. [8], where
distances from jamming are defined by the averaged overlaps. In this study, we prepare
diﬀerent 50 packings for each distance from jamming, ∆ϕ ≡ ϕ−ϕJ , by changing the initial
random configurations of the particles, where ϕ and ϕJ are defined as the area fraction of
the particles and the jamming point [3], respectively.
We apply a simple shear deformation to the static packings by replacing every particle
positions, (xi , yi ), with (xi + ∆γyi , yi ) (i = 1, . . . , N ), where ∆γ is an infinitesimal strain
increment and the Lees-Edwards boundary conditions are used.
3

A Master Equation

Macroscopic mechanical responses of granular materials to simple shear deformations
are described by constitutive equations [1]. In this section, we show that constitutive
equations can be described by microscopic mechanical responses of force-chain networks
through a master equation for the PDFs of forces. In two-dimensional soft particle packings, a microscopic expression of stress tensor is given by σpq = S −1 ⟨fp dq ⟩ (p, q = x, y),where S, fp , and dq are the system area, interparticle force, and relative position (orbranch vector) between two particles in contact, respectively. On the right-hand-side,
the parentheses, ⟨. . . ⟩, represent a statistical average over the particles in contacts. In
our linear frictionless model, the force is proportional to the overlap and parallel tothe relative position. Therefore, introducing a unit vector along the relative positionas (nx , ny ) = (cos θ, sin θ) with an angle measured from the x-axis (Fig. 1(a)), we canrewrite the force and relative position as fp = Kξnp and dq = d(ξ)nq , respectively, whered(ξ) = s − ξ is the interparticle distance defined as the diﬀerence between the sum ofradii, s, and overlap. Then, the stress tensor is rewritten as
σpq =K⟨ξd(ξ)np nq ⟩ ,S(1)

where the stress tensor is represented by a set of overlap and angle, i.e. r ≡ (ξ, θ), which
we call microscopic states of soft particle packings.
2

Kuniyasu Saitoh and Fumiko Ogushi

Statistics of the microscopic states are governed by their PDFs, Pγ (r), where the subscript, γ, represents the amount of shear strain applied to the system. To introduce thePDFs, we employ the Delaunay triangulations to static packings, where not only the particles in contacts, but also the nearest neighbors without contacts, i.e. the particles invirtual contacts, are connected by the Delaunay edges. Then, we generalize “overlaps”
as ξ ≡ s − D with the Delaunay edge length, D, where the overlaps between particles invirtual contacts are defined as negative values (s < D). Because the total number of theDelaunay edges is conserved during deformations, the PDFs are normalized as∫∫Pγ (r)dr = 1 .
(2)
When the system is deformed by a small shear strain, ∆γ, the microscopic states willchange to new states, r′ ≡ (ξ ′ , θ′ ), as shown in Fig. 1(b), where the new PDF is given byPγ+∆γ (r′ ). Assuming that the transitions between microscopic states (from r to r′ ) aredescribed by the Markov processes, we connect the PDFs before and after simple shear
deformation through the Chapman-Kolmogorov equation [10],∫∫′Pγ+∆γ (r ) =W (r′ |r)Pγ (r)dr ,(3)
where W (r′ |r) is the CPD of the microscopic states, r′ , which were r before the simple
shear deformation. By definition, the CPD is normalized as∫∫
W (r′ |r)dr′ = 1 .
(4)
From Eqs. (3) and (4), a master equation for the PDFs is readily obtained as∫∫∂′Pγ (r ) ={T (r′ |r)Pγ (r) − T (r|r′ )Pγ (r′ )} dr ,∂γ

(5)where we introduced the transition rate as T (r′ |r) ≡ lim∆γ→0 W (r′ |r)/∆γ [10]. In themaster equation (5), the first and second terms in the integral on the right-hand-side
represent the gain and loss of new states, r′ , respectively. Therefore, the transition rates
or CPDs fully determine the statistics of microscopic responses to shear deformations.
If we take the statistical average, ⟨. . . ⟩ in Eq. (1), by the PDF, the shear stress definedas τ = (σxy + σyx )/2 is given by∫∫Kτ={ξd(ξ) sin 2θ}Pγ (ξ, θ)dξdθ ,(6)2S
where we used nx ny = cos θ sin θ = (sin 2θ)/2 and restricted the overlaps to positive values(ξ > 0). Then, the shear modulus, G = ∂τ /∂γ, is given by
{
}
∫∫
∂
K
{ξd(ξ) sin 2θ}Pγ (ξ, θ) dξdθ ,
(7)
G=
2S
∂γ
which is coupled with the master equation (5) through the derivative of the PDF.
3

Kuniyasu Saitoh and Fumiko Ogushi

4

Results

4.1 Microscopic responses
To determine the transition rates in the master equation, we first study aﬃne and nonaﬃne responses of microscopic states, r = (ξ, θ). When we apply simple shear deformation
to the system, the microscopic state will change to raﬃne ≡ (ξaﬃne , θaﬃne ), where eachcomponent is given by the initial state asξaﬃne ≃ ξ −∆γd(ξ) sin 2θ ,2

θaﬃne ≃ θ +∆γ
cos 2θ ,2(8)

respectively 1 . The second term on the right-hand-side of ξaﬃne represents an anisotropic
response of overlaps, where its magnitude becomes maximum in compressive and decompressive directions, θ = ∓π/4, respectively. Note that Eq. (8) is correct for both contacts
and virtual contacts. In addition, if we introduce change rates of microscopic states asṙaﬃne ≡ (raﬃne − r)/∆γ, each component is respectively given byd(ξ)ξ˙aﬃne = −sin 2θ ,2

θ̇aﬃne =1
cos 2θ .
2(9)

In soft particle packings, however, particles are randomly arranged and the force balance between the particles in contacts is broken by simple shear deformation so thatparticles move around and the system relaxes to a new static state, where the microscopic
states change to r′ = (ξ ′ , θ′ ). Clearly, the microscopic states after relaxation, r′ , as well˙ θ̇), are expected to be diﬀerent from those foras their change rates, ṙ ≡ (r′ − r)/∆γ ≡ (ξ,aﬃne response, Eqs. (8) and (9). Figures 1(c) and (d) display the angular dependence ofchange rates, where the aﬃne responses (blue dots) show the deterministic evolution according to Eq. (9), while the non-aﬃne responses (red dots) fluctuate around mean values.
˙ seems to be quite diﬀerent fromIn these figures, the mean change rate of overlaps, ⟨ξ⟩,˙ξaﬃne , while that of angles is almost the same with θ̇aﬃne , i.e. ⟨θ′ ⟩ ≃ θaﬃne . Therefore, we
can write the non-aﬃne microscopic responses as the sum of their mean and fluctuations,i.e.
ξ ′ = f (ξ, θ) + ψ , θ′ = θaﬃne + ζ ,(10)respectively, where the mean overlap, ⟨ξ⟩ ≡ f (ξ, θ), depends on the initial states andwe have √introduced random√ variables as ψ and ζ of which standard deviations are given
by v ≡ ⟨ψ 2 ⟩ and w ≡ ⟨ζ 2 ⟩, respectively. As can be seen in Figs. 1(c) and (d), thestandard deviations, v and w, are independent of the initial angle, i.e. they are isotropic
in space.
1

Eq. (8) is the first-order approximations of an infinitesimal shear strain, O(∆γ), where the frame
rotation by simple shear deformation, φ = −∆γ/2, has already been subtracted from θaﬃne and θ′ .

4

Kuniyasu Saitoh and Fumiko Ogushi

(a)

(c)

(d)

2

4

1

2

0

0

-1

-2

(b)

-2
-π/2

-π/4

0

π/4

π/2

-4
-π/2

-π/4

0

π/4

π/2

Figure 1: (Color online) Sketches of microscopic states (a) before and (b) after simple shear deformation,r = (ξ, θ) and r′ = (ξ ′ , θ′ ), respectively, where the solid lines connect centers of particles and the dotted˙lines are parallel to the x-axis. (c) Change rates of overlaps scaled by interparticle distances, ξ/d,and (d)those of angles, θ̇, plotted against the initial angle, θ, where the blue and red dots are the results of aﬃne
and non-aﬃne responses, respectively. The dotted lines in (c) and (d) are − sin(2θ)/2 and cos(2θ)/2,
respectively. Here, we apply ∆γ = 10−6 to the system with ∆ϕ = 10−2 (i.e. α ≡ ∆γ/∆ϕ = 10−4 ).

4.2 Mean and fluctuations
The mean overlap, f (ξ, θ), can be determined by scattered plots of overlaps before andafter simple shear deformation. Figure 2 displays the scattered plots, where we plot theaﬃne and non-aﬃne overlaps after shear, ξaﬃne (blue dots) and ξ ′ (red dots), against theinitial overlap, ξ. Here, we fixed the initial angle, θ, to compressive (Figs. 2(a) and (b)) or
decompressive direction (Figs. 2(c) and (d)). The diﬀerence between aﬃne and non-aﬃneresponses depends on both the strain step, ∆γ, and distance from jamming, ∆ϕ, where thediﬀerence increases with their ratio, α ≡ ∆γ/∆ϕ. Because overlaps are defined not only
between contacts, but also between virtual contacts, there are four kinds of transitions, i.e.
contact-to-contact (ξ, ξ ′ > 0), contact-to-virtual (ξ > 0, ξ ′ < 0), virtual-to-contact (ξ < 0,
ξ ′ > 0), and virtual-to-virtual (ξ, ξ ′ < 0) [8]. In this paper, we restrict our analyses
to the first case, contact-to-contact, where both overlaps before and after simple shear
deformation remain positive. Then, we introduce the mean overlap for contact-to-contact
as a linear fitting function for ξ ′ (> 0) asf (ξ, θ) = {1 + a(θ)} ξ + b(θ) .

(11)We also introduce a standard deviation of ξ ′ around f (ξ, θ) as v such that systematicdeviations from aﬃne responses are quantified by the excess slope, a(θ), oﬀset, b(θ),
and fluctuation, v, where the aﬃne approximation gives a(θ) = (∆γ/2) sin 2θ, b(θ) =−(s∆γ/2) sin 2θ, and v = 0.
Figures 2(e) and (f) display our numerical results of the excess slopes and oﬀsets, where
their angular dependence are well captured by sinuous functions,a(θ) = −A sin 2θ ,b(θ) = −B sin 2θ ,5

(12)

Kuniyasu Saitoh and Fumiko Ogushi

6

2

6

(a)

(b)

1x10-2(1x10-3)
3x10-3(3x10-3)
1x10-3 (1x10-2)-4-2
3x10 (3x10 )
1x10-4 (1x10-1)

(e)

3

3

1

0

0

0
2

1

-3

-3

-1

0

-1

-2
-π/2

-6

-6
-6

-3

0

3

6

6

-6

-3

0

3

6

6

-π/4

-2
-π/2

0

π/4

π/2

-π/4

0

4

(c)

(d)
3

2

0

0

0

π/2

1x10-2(1x10-3)
3x10-3(3x10-3)
1x10-3(1x10-2)-2
3x10-4(3x10)
1x10-4 (1x10-1)

(f)

3

π/4

2

1

-3

-3

-2

0

-1

-6

-2
-π/2

-6
-6

-3

0

3

6

-6

-3

0

3

6

-4
-π/2

-π/4

0

-π/4

π/4

π/2

0

π/4

π/2

Figure 2: (Color online) Scatter plots of ξ ′ against ξ scaled by the averaged overlap in the system, where
we fix the initial angle to θ = −π/4 ((a) and (b)) or π/4 ((c) and (d)), and change the ratio, α, as shownin the legends. The angular dependence of (e) excess slopes, a(θ), and (f) oﬀsets, b(θ), where the lines are
fitting functions, (e) −A sin 2θ and (f) −B sin 2θ, respectively. The insets show the data scaled by each
amplitude, i.e. a(θ)/A and b(θ)/B, where the dotted lines represent − sin 2θ. Here, we apply ∆γ = 10−5to the system with diﬀerent ∆ϕ (or the ratio, α ≡ ∆γ/∆ϕ) as shown in the legends.

respectively. On the right-hand-side of Eq. (12), the amplitudes, A and B, represent thestrength of anisotropic responses of soft particle packings. We have confirmed that theamplitudes, A and B, and fluctuation, v, are increasing functions of the ratio 2 , α. From
Eq. (12), the mean overlap, Eq. (11), is rewritten asf (ξ, θ) = ξ − (Aξ + B) sin 2θ ,(13)

where the sinuous dependence is the same with that of aﬃne response, Eq. (8).
4.3 Conditional probability distributions
The microscopic aﬃne responses are described by the deterministic equation (8) sothat the CPDs are given by a delta-function, W (r′ |r) = δ(r − raﬃne ). However, as shownin Figs. 1 and 2, the microscopic non-aﬃne responses fluctuate around their mean values
such that the CPDs must have finite widths. In this section, we numerically determine
the CPDs by MD simulations.
2

Data are not shown. See Ref. [8] for the case of isotropic compressions.

6

Kuniyasu Saitoh and Fumiko Ogushi

(a)

3

10

(b)

-3

(c) 103

3

10

0
|∆γ|
2|∆γ|
3|∆γ|
4|∆γ|

1x10-3
3x10
-2
1x10

10

2

10

2

102

1

10

1

101

10

0

10
-0.04

0

-0.02

0

0.02

0.04

10
-0.008-0.004

0

0.004

0.008

100
-0.008

-π/3
-π/4
-π/6
0
π/6
π/4
π/3

-0.004

0

0.004

0.008

Figure 3: (Color online) CPDs for contact-to-contact, W (r′ |r), plotted against the distance from themean overlap, Ξ ≡ ξ ′ − f (ξ, θ), which is scaled by the averaged overlap in the system, where we fixed theinitial overlap to ξ = 1.2. Diﬀerent symbols represent diﬀerent (a) α, (b) ∆θ, and (c) θ, as shown in thelegends, where we used (a) θ = −π/4 and ∆θ = 0, (b) α = 10−3 and θ = −π/4, and (c) α = 10−3 and∆θ = 0, respectively. The lines are the Gaussian fitting function, Eq. (14).

Figure 3 displays the CPDs for contact-to-contact, W (r′ |r), plotted against the distancefrom the mean overlap, Ξ ≡ ξ ′ − f (ξ, θ). In Fig. 3(a), we fix the initial angle and angle
diﬀerence to θ = −π/4 and ∆θ ≡ θ′ − θaﬃne = 0, respectively, while we change the ratio
from α = 10−3 to 10−2 . In Fig. 3(b), the initial angle and ratio are fixed to θ = −π/4 andα = 10−3 , respectively, while the angle diﬀerence changes from |∆θ| = 0 to 5∆γ (i.e. from|θ̇| = 0 to 5). In Fig. 3(c), the ratio and angle diﬀerence are given by α = 10−3 and ∆θ = 0,
respectively, but we change the initial angle from θ = −π/3 to π/3. From these results,
the CPDs are symmetric around the mean value and change their shapes depending on
the ratio, α, and angle diﬀerence, ∆θ, while they are independent of the initial angle,
i.e. the CPDs are isotropic implying that fluctuations of forces are equiprobable in all
directions. Therefore, we can decompose the CPDs as
W (r′ |r) = C(θ′ |θ)Z(ξ ′ |r) ,(14)

where C(θ′ |θ) is another CPD of θ′ which was θ before simple shear deformation. Note that
C(θ′ |θ) depends only on the angle diﬀerence, ∆θ, but not on the initial angle, θ, i.e. C(θ′ |θ)is also isotropic, implying that the probability of angle fluctuations is also equiprobablein all directions. The θ′ -independent function is given by a Gaussian distribution,
1
2
2
Z(ξ ′ |r) = √e−Ξ /2v ,
2πv

(15)

where the lines in Fig. 3 represent Eq. (14) with the Gaussian distribution, Eq. (15).
5Summary

In this study, we have investigated microscopic mechanical changes of granular materials under simple shear deformations by the master equation and MD simulations. We havefound that, though microscopic aﬃne responses are deterministic, microscopic non-aﬃne
7

Kuniyasu Saitoh and Fumiko Ogushi

responses fluctuate around mean values as a result of complicated non-aﬃne displacements
during relaxation. Such the stochastic evolution or restructuring of force-chain networks
is well captured by the mean overlap, Eq. (13), and CPD, Eq. (14), where anisotropic
mechanical responses of granular materials are well described by the sinuous functions
for the excess slope and oﬀset, Eq. (12), and the CPD for contact-to-contact is foundto be the isotropic Gaussian distribution function, Eq. (15). Therefore, the fluctuationsof overlaps (forces) and angles are equiprobable in all directions even though mechanicalresponses of granular materials are anisotropic in average.
Acknowledgment
This work was financially supported by World Premier International Research Center
Initiative (WPI), MEXT, Japan. A part of numerical 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] M. van Hecke, J. Phys.: Condens. Matter (2010) 22, 033101.
[4] S. Henkes and B. Chakraborty, Phys. Rev. E (2009) 79, 061301.
[5] J. H. Snoeijer, T. J. H. Vlugt, M. van Hecke, and W. van Saarloos, Phys. Rev. Lett.
(2004) 92, 054302.
[6] E. I. Corwin, H. M. Jaeger, and S. R. Nagel, Nature (2005) 435, 1075.
[7] L. E. Silbert, G. S. Grest, and J. W. Landry, Phys. Rev. E (2002) 66, 061303.
[8] K. Saitoh, V. Magnanimo, and S. Luding, Soft Matter (2015) 11, 1253.
[9] S. Luding, J. Phys.: Condens. Matter (2005) 17, S2623.
[10] N. G. van Kampen, Stochastic Processes in Physics and Chemistry, 3rd edition(2007). Elsevier B. V. Amsterdam, The Netherlands.

8