Influence of polydispersity on micromechanics of granular materials

# Influence of polydispersity on micromechanics of granular materials

M. Reza Shaebani Department of Theoretical Physics, University of Duisburg-Essen, 47048 Duisburg, Germany    Mahyar Madadi Department of Applied Mathematics, The Research School of Physics and Engineering, The Australian National University, Canberra 0200, Australia Department of Exploration Geophysics, Curtin University, P.O. Box U1987, Perth, WA 6845, Australia    Stefan Luding Multi Scale Mechanics (MSM), CTW, UTwente, P.O. Box 217, 7500 AE Enschede, Netherlands    Dietrich E. Wolf Department of Theoretical Physics, University of Duisburg-Essen, 47048 Duisburg, Germany
July 5, 2019
###### Abstract

We study the effect of polydispersity on the macroscopic physical properties of granular packings in two and three dimensions. A mean-field approach is developed to approximate the macroscale quantities as functions of the microscopic ones. We show that the trace of the fabric and stress tensors are proportional to the mean packing properties (e.g. packing fraction, average coordination number, and average normal force) and dimensionless correction factors, which depend only on the moments of the particle-size distribution. Similar results are obtained for the elements of the stiffness tensor of isotropic packings in the linear affine response regime. Our theoretical predictions are in good agreement with the simulation results.

###### pacs:
45.70.-n, 45.70.Cc, 83.80.Fg

## I Introduction

The physics of granular media has received a lot of attention because of its scientific challenges and industrial relevance. The structural and dynamical properties of granular materials differ from those of ordinary solids, liquids, or gases due to nonlinearity and disorder deGennes99 (); Jaeger96 (); Hinrichsen04 (). On the microscopic level, a static assembly of grains consists of particles which interact with their neighbors in order to prevent interpenetration. In spite of the uniform density of granular packings, the resulting contact and force networks between particles are highly inhomogeneous Radjai96 (); Makse00 (); Ostojic06 (), leading to many intriguing phenomena in these systems. Describing the behavior via micromechanical approaches, in which the discrete nature of the system is taken into account, is thus commonly preferred to continuum-mechanical approaches where some heuristic assumptions have to be made in order to construct the constitutive equations for macroscopic fields. One can then express the macroscopic physical quantities in terms of the microscale ones. For example, thermal and electrical conductivities are related to the trace of the fabric tensor, a micro-geometrical probability of the orientations of contacts. While the relationship between macroscopic and microscopic properties of granular media has been studied widely deGennes99 (); Hinrichsen04 (); Luding04 (), the question remains as to what extent the macroscale quantities are sensitive to the micro-scale details, and how large is the error introduced in the calculation of the “observable quantities” by taking into account only the average packing properties.

Granular materials in nature and industry consist of particles with the common property of polydispersity. It is known that size polydispersity affects the mechanical behavior of granular systems (e.g. shear strength) Wackenhut05 (); Voivret09 () as well as their space-filling properties (e.g. packing fraction) Herrmann03 (); Voivret07 (), which are crucial in many chemical processes like absorption, filtering, etc. Polydispersity in most studies, so far, has been restricted to narrow size distributions mainly to prevent long-range structural order; however, there are a few studies where broader ranges of particle-size distribution are investigated Voivret09 (); Voivret07 (); Dodds02 (); Ogarko12 (). In this paper, we address the question of how deviation from the monodisperse case influences the macroscopic properties of granular assemblies.

We consider a special case of spherical particles [or disks in two dimensions (2D)] allowing for analytical calculations. The main goal is to develop a mean-field approach to calculate the desired microscopic quantities, such as the trace of the fabric and stress tensors, and the elements of the stiffness tensor in two- and three-dimensional polydisperse granular systems. These quantities are directly connected to macroscopic quantities such as thermal and electrical conductivities, isotropic pressure, and bulk and shear moduli. A similar analytical approach has been already used in Ref. Madadi04 () to calculate the trace of the fabric tensor in 2D packings, where it turned out that the trace of fabric is factorized into three contributions: (i) the volume fraction, (ii) the mean coordination number, and (iii) a dimensionless correction factor which only depends on the particle-size distribution. Using a similar approach, here we investigate also the stress and stiffness tensors and extend the method to 3D cases. In order to compare the analytical results with numerical simulations, we first construct static packings of grains using contact dynamics simulations Moreau94 (); Jean99 (); Brendel04 (). The initial dilute systems of rigid particles are compressed by imposing a confining pressure to get the final static homogeneous packings Shaebani09 (). Comparisons are then made between the results of our mean-field model and the exact values obtained from the numerical simulations.

This work is organized in the following manner: The fabric tensor of a polydisperse assembly of spherical particles is investigated in Sec. II, and a mean-field approach is introduced to calculate the trace of fabric. We present the analytical results for the calculation of the stress tensor in Sec. III, and the same approach is used in Sec. IV to investigate the stiffness tensor elements in frictionless isotropic packings. In Sec. V, the analytical calculations are compared to numerical simulations of corresponding packings of polydisperse particles. Finally, we discuss and conclude the results in Sec. VI. Detailed calculations for two-dimensional packings of disks are presented in the Appendix.

## Ii Fabric tensor

### ii.1 Single-particle case

Various definitions of the fabric tensor have been used in the literature to describe the spatial arrangement of the particles in a granular assembly Goddard86 (); Chang88 (); Rothenburg81 (). The fabric tensor of the second order for one particle is defined as Mehrabadi82 (); Latzel00 (); Madadi04 ()

 hpαβ=Cp∑c=1lpcα|→lpc|lpcβ|→lpc|, (1)

where is the number of contacts of particle , and is the component of the branch vector , connecting the center of particle to its contact . In the case of spherical particles, the unit branch vector and the unit normal vector at contact are identical. The trace of the single-particle fabric tensor in a -dimensional system is

 hpαα=Cp∑c=1D∑α=1lpcα|→lpc|lpcα|→lpc|=Cp, (2)

i.e. the number of contacts of particle .

### ii.2 Many-particle case

The average fabric tensor enables us to describe the global contact network in a given volume . Assuming that the contribution of particle (lying inside ) to the average fabric tensor is proportional to its volume , we obtain

 ⟨hαβ⟩V=1VN∑p=1Vphpαβ, (3)

where the sum runs over all particles lying inside , and denotes the volume weighted average. Using Eq. (2) to calculate the trace of the average fabric tensor, we get

 ⟨hαα⟩V=1VN∑p=1VpCp, (4)

which can be interpreted as the contact number density. Alternative possibilities, e.g. using the volume of the polygon that contains the particle (obtained e.g. via Voronoi tessellation), or introducing constant prefactors or slightly different volume contributions are not discussed here (see Refs. Goddard86 (); Chang88 (); Cowin85 (); Duran10 () for more details). In a monodisperse packing, Eq. (4) for identical particles is reduced to , where is the packing fraction (), and is the average coordination number (). We note that only “real” contacts contribute to the calculation of , and geometrical neighbors without a permanent physical contact, which do not contribute in the fabric and force carrying structures, are not considered here.

### ii.3 Polydispersity

For an accurate evaluation of the trace of the average fabric tensor in a polydisperse granular packing, one should take into account the contributions from all particles. However, if the distribution function of particle radii is known, can be approximated as a function of the moments of the size distribution. We assume a polydisperse distribution of particle radii with probability to find the radius between and , and with . The continuum limit of Eq. (4) is then given by

 ⟨hαα⟩V=NV∫∞0V(a)C(a)f(a)da. (5)

Here, is the average coordination number of particles with radius . We evaluate using a mean-field approach similar to the one proposed in Ouchiyama81 () and used already in Madadi04 () to study the trace of the fabric tensor. In the following, we concentrate on the case of spherical particles in three-dimensional systems (see the detailed calculations for two-dimensional packings of disks in the Appendix). Let us suppose that each particle in the polydisperse granular medium is surrounded by identical particles of average radius (see Fig. 1), where . The surface of a reference particle of radius is then shielded by its neighboring particles of radius . The space angle covered by a neighboring particle on the reference particle in a three-dimensional packing of spheres is

 Ω(a)=2π(1−√(a+⟨a⟩)2−⟨a⟩2a+⟨a⟩). (6)

The total fraction of shielded surface, also called linear compacity, is obtained as

 cs(a)=14πa2C(a)∑i=1Ω(a)a2=Ω(a)C(a)/4π. (7)

Now, another basic assumption is that the total fraction of shielded surface is independent of the particle radius . As a result, the expected mean coordination number becomes

 z=∫∞0C(a)f(a)da=4πcsq0, (8)

with . Using Eqs. (7) and (8) one finds

 C(a)=zq0Ω(a). (9)

The trace of the fabric tensor for a polydisperse packing is then obtained by substitution of Eq. (9) in Eq. (5),

 ⟨hαα⟩V=ϕzg1, (10)

where the correction factor is defined as

 g1=∫∞0V(a)f(a)Ω(a)daq0∫∞0V(a)f(a)da=⟨a3⟩g⟨a3⟩ (11)

Here, and denote the th moments of the size distribution and the modified distribution normalized by , respectively. We note that depends only on the size distribution function .

### ii.4 Narrow size distributions

By introducing , which ranges between and depending on the choice of , Eq. (6) can be written as

 Ω(a)=2π(1−√ϵ2+4ϵ+32+ϵ). (12)

Indeed, quantifies the deviation from the mean particle size [e.g.  equals zero in the monodisperse case]. Hence, for narrow size distributions, we approximate by Taylor expansion around (corresponding to Taylor expansion around ). By Taylor expansion to the second order in , one obtains

 1Ω(a)≃A1+B1ϵ+C1ϵ2, (13)

with , , and . The first-order approximation deviates significantly from the exact value (see Fig. 2). However, the second-order expansion provides a good approximation with less than error in the range (or ).

Therefore, the correction factor [Eq. (11)] for narrow size distributions becomes

 (14)

Equation (14) should account for arbitrarily shaped size distributions as long as they are not too wide. Note the different nomenclature in Ref. Goncu10 (), where the above equation is introduced with different abbreviations and coefficients.

## Iii Stress tensor

### iii.1 Single-particle case

The micromechanical expressions for the components of the stress tensor of a single particle in a static granular assembly are Christoffersen81 (); Latzel00 ()

 σpαβ=1VpCp∑c=1lpcαFpcβ, (15)

where is the force exerted on particle by its neighboring particle at contact .

One could assume in a crude approximation that the force at contact is equal to in a three-dimensional system, where , , and are the average normal and tangential contact forces around the particle , and , , and are the normal and tangential unit vectors at contact , respectively. Then the force-averaged stress tensor becomes

 (16)

For a spherical grain, we project the contact unit vectors (, , ) onto an arbitrary Cartesian coordinate system [Figs. 3(a) and 3(b)], and write the force-averaged stress tensor of a single particle as

 ˜σp=apVpCp∑c=1[¯Fpn⎛⎜⎝W2002W2011W1101W2011W2020W1110W1101W1110W0200⎞⎟⎠ +¯Fpt1⎛⎜⎝W1102W1111−W2001W1111W1120−W2010W0201W0210−W1100⎞⎟⎠ +¯Fpt2⎛⎜⎝−W1011W10020−W1020W10110−W0110W01010⎞⎟⎠], (17)

where the function is defined as

 Wmnkl=sinm(θc)cosn(θc)sink(φc)cosl(φc), (18)

with and . Using Eq. (16), the trace of the stress tensor becomes

 ˜σpαα=apVpCp∑c=13∑α=1(¯Fpnnpcαnpcα+¯Fpt1npcαtpc1α+¯Fpt2npcαtpc2α) =apVpCp∑c=1(¯Fpn|^npc|2+¯Fpt1^npc⋅^tpc1+¯Fpt2^npc⋅^tpc2)=apVp¯FpnCp.

Equation (LABEL:stress-single-trace) remains valid also in the 2D case (see Appendix). As expected for isotropic packings, the trace of the stress tensor and therefore the isotropic pressure () do not depend on the tangential forces.

### iii.2 Many-particle case

In the many-particle case, the average stress tensor in a given volume is defined as Latzel00 ()

 ⟨σαβ⟩V=1VN∑p=1Vpσpαβ=1VN∑p=1Cp∑c=1lpcαFpcβ, (20)

where the sum runs over all particles lying inside . Using Eq. (LABEL:stress-single-trace) to calculate the trace of the average stress tensor, we get

 ⟨˜σαα⟩V=1VN∑p=1Vp˜σpαα=1VN∑p=1ap¯FpnCp. (21)

### iii.3 Polydispersity

Now we assume a polydisperse distribution of particle radii with probability to find the radius between and , and with . Assuming that the average contact force exerted on a particle depends only on its radius , the continuous limit of Eq. (21) in a mean-field approximation is given by

 ⟨˜σαα⟩V=NV∫∞0a¯Fn(a)C(a)f(a)da. (22)

In Eq. (22), it is supposed that all particles of size have a certain mean coordination number and a certain mean normal force . Indeed, particles of the same size may have different coordination number and normal contact forces, however, the main goal here is to propose a method to calculate macroscopic quantities without taking into account all the microscopic details of the system. We use the mean-field approach introduced in Sec. II.3 to evaluate . By substitution of Eq. (9) in Eq. (22) we get

 ⟨˜σαα⟩V=NVz∫∞0a¯Fn(a)f(a)Ω(a)daq0 (23) =ϕz∫∞0a¯Fn(a)f(a)Ω(a)daq0∫∞0V(a)f(a)da.

According to the mean-field approach used in Sec. II.3, increases with increasing radius . Now, let us assume that the average normal force also increases with , so that the ratio remains roughly constant Madadi05 (). We calculate this ratio for the average-sized particles in the following:

 ¯Fn(a)C(a)=¯Fn(a)zq0Ω(a)≃¯Fn(⟨a⟩)zq0Ω(⟨a⟩)=q0¯Fn(⟨a⟩)Ω(⟨a⟩)z, (24)

therefore

 ¯Fn(a)=Ω(⟨a⟩)¯Fn(⟨a⟩)Ω(a). (25)

By substitution of Eq. (25) in Eq. (23), we obtain

 ⟨˜σαα⟩V=3ϕz¯Fn(⟨a⟩)g24π⟨a2⟩ (26)

with

 g2=(2−√3)π⟨a2⟩∫∞0af(a)Ω2(a)daq0⟨a3⟩ (27)

### iii.4 Narrow size distributions

In the limit of narrow size distributions, we approximate by Taylor expansion around (similar to Sec. II.4):

 1Ω2(a)≃A2+B2ϵ+C2ϵ2, (28)

with , , and . By substitution of Eq. (28) in Eq. (27) we obtain the correction factor for arbitrary narrow distributions,

 (29)

## Iv Stiffness tensor

The linear response of a material to “weak” external perturbations is described by a fourth rank tensor, which is called the elastic or stiffness tensor Kruyt96 (); Luding04 (). This tensor has and elements in three- and two-dimensional systems, respectively, but they are not all independent. Symmetry considerations reduce the number of independent elements. For example, the elastic behavior of isotropic materials can be described by only two independent parameters, usually represented by Lamé coefficients and . In this section, the stiffness tensor of a homogeneous and isotropic assembly of polydisperse particles is investigated (for the case of an anisotropic monodisperse system, see, e.g. Shaebani11 (); Mouraille06 ()).

The stiffness tensor for a spherical particle, where affine deformation is assumed, is defined as Bathurst88 (); Luding04 ()

 Cpα,β,γ,η=2a2pVpCp∑c=1(knnpcαnpcβnpcγnpcη+ktnpcαtpcβnpcγtpcη), (30)

where is the unit vector parallel to the tangential component of the contact force [see Fig. 3(c)]. The volume weighted average of is then given by

 ⟨Cα,β,γ,η⟩V=1VN∑p=1VpCpα,β,γ,η= (31) 1VN∑p=12a2pCp∑c=1(knnpcαnpcβnpcγnpcη+ktnpcαtpcβnpcγtpcη).

Note that the stiffness tensor is basically determined by the packing geometry. For ease of calculation, we consider only frictionless packings, i.e.  is set to zero hereafter. Using the microscopic information of the contact orientations, one can accurately calculate the elements of via Eq. (31). Next, the Lamé constants and can be deduced from the stiffness tensor, e.g. as and or, more generally, as and where

 ⟨¯¯¯¯¯¯¯Ciijj⟩V=1D(D−1)D∑i≠j⟨Ciijj⟩V,⟨¯¯¯¯¯¯¯Ciiii⟩V=1DD∑i⟨Ciiii⟩V, (32)

and is the dimension of the system. The macroscopic physical quantities of interest are the bulk modulus and the shear modulus , which can be deduced from the Lamé coefficients in isotropic materials as

 G/kn=μ/kn=⟨¯¯¯¯¯¯¯Ciiii⟩V−⟨¯¯¯¯¯¯¯Ciijj⟩V2kn, (33)

and

 K/kn=(λ+2Dμ)/kn=⟨¯¯¯¯¯¯¯Ciiii⟩V+(D−1)⟨¯Ciijj⟩VDkn. (34)

Now, assuming a polydisperse probability distribution of particle radii , Eq. (31) for can be written as

 ⟨Cα,β,γ,η⟩V=NknV∫∞02a2(C(a)∑c=1ncαncβncγncη)f(a)da. (35)

Since the packings are supposed to be isotropic and homogeneous, we assume that grains are scattered homogeneously around the reference particle. Therefore, the summation over neighbors can be approximated by the following integration in three dimensions (for the 2D case, see Appendix):

 C(a)∑c=1Q(θc,φc)=C(a)4π∫π0dθsin(θ)∫2π0dφQ(θ,φ). (36)

We present the reduced form of the fourth rank tensor by mapping , i.e. , , , , , and . Using Eqs. (9), (35), and (36), one obtains

 ⟨C⟩V=Nknz2πVq0∫∞0daa2Ω(a)f(a)× ∫π0dθsin(θ)∫2π0dφ⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝W4004W4022W2202W4013W3103W3112W4040W2220W4031W3121W3130W0400W2211W1301W1310W4022W3112W3121W2202W2211W2220⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠,

where elements were defined in Eq. (18). After integration on and , the volume weighted average of the stiffness tensor for an isotropic polydisperse packing becomes

 ⟨C⟩V=ϕzkng310π⟨a⟩⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝311000310003000100101⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠. (38)

The correction factor is defined as

 g3=⟨a⟩⟨a2⟩g⟨a3⟩, (39)

and for narrow size distributions, one obtains

 g3≃(A1−B1+C1)⟨a⟩⟨a2⟩⟨a3⟩+(B1−2C1)+C1⟨a4⟩⟨a⟩⟨a3⟩(A1−C1)+C1⟨a2⟩⟨a⟩2, (40)

with the same coefficients as defined after Eqs. (13) and (28). To summarize this section, the Lamé constants for frictionless isotropic packings are

 μ=λ=(knϕzg3)/(10π⟨a⟩), (41)

and the shear and bulk moduli are

 G/kn=(ϕzg3)/(10π⟨a⟩), (42)

and

 K/kn=(ϕzg3)/(6π⟨a⟩). (43)

Notably, in three-dimensional frictionless isotropic packings, independent of their size distribution and average packing properties.

## V Simulation results

To verify the theoretical predictions of the previous sections, we carry out numerical simulations with the help of the contact dynamics (CD) algorithm Moreau94 (); Jean99 (); Brendel04 (). We first construct 2D and 3D static homogeneous packings in zero gravity by compressing the initial dilute configuration of particles [Fig. 4 (left)]. Periodic boundary conditions are imposed in all directions to avoid the side effects of lateral walls. The compaction is achieved by imposing a constant external pressure and letting the size of the system evolve in time Andersen80 (). As the volume of the system decreases, after a while particles touch each other and build an inner pressure , which resists and eventually compensates , so that finally equals . Particles prevent further compaction, and a static homogeneous configuration is reached [Fig. 4 (right)]. The full description of the packing generation method can be found in Shaebani09 (). In order to illustrate the validity range of our assumptions, we generate three types of polydisperse packings with uniform particle-size distributions but with different widths (see Table LABEL:Table1). We denote the samples SMP1, SMP2, and SMP3, respectively, by full circles, open squares, and full triangles throughout this section. To investigate the effect of friction, we construct a new packing for each value of the particle-particle friction coefficient . In particular, the results corresponding to , , and are hereafter denoted by green, blue and red colors. The number of grains contained by packings are and in 2D and 3D cases, respectively.

For comparison with the theory, we first test the validity of assumptions made in Sec. II.3. The linear compacity is displayed in Fig. 5 for the static configurations of particles obtained from the isotropic compression simulations. For each particle , the surface angle covered by its neighboring particle at contact is calculated, and the linear compacity of particle is obtained as or for two- or three-dimensional packings, respectively. Next, we divide the range of possible values of the particle radius into bins. Each data point in Fig. 5 corresponds to the mean value of , averaged over all particles in the same bin. The contribution of the rattler particles, which transmit no force, is excluded. For moderate widths of size distributions (SMP1), is approximately constant in for a given packing (we note that the fluctuations of around its mean value in a given packing originate from the finite size of the samples). However, is remarkably above the average value for small particle sizes in wider distributions (SMP2 and SMP3). This is a common property of our highly polydisperse packings (with uniform size distribution) that the fraction of shielded surface is larger than the average for small particles if rattlers are excluded (see Ogarko12 () for uniform volume distributions). A similar behavior has been observed in discrete element method simulations of soft particles Goncu10 (). There, it is also shown that if rattlers are included in the statistics, the small particles on average are less covered than the larger ones. However, the deviation of small particles from the average decreases as the volume fraction of the packing increases by incremental compression.

Another point is that depends strongly on the dimension of the system and the friction coefficient. Increasing the friction stabilizes the system in a less dense state and decreases the connectivity of the contact network Kadau03 (); Shaebani09b (). Therefore, we expect lower values of and when increasing , as confirmed by the data.

In Figs. 6 and 7, the coordination number is shown as a function of for the same set of systems as in Fig. 5. For comparison, we also plot from Eq. (9). Here, the average coordination number of the packing is taken from the simulation results, is provided by Eq. (6) or (44), and the size distribution of each packing after the compaction process is used to calculate . The mean-field approach of Sec. II.3 qualitatively fits well to the data, however, the slopes of the curves are slightly greater than the corresponding slopes of the best-fit curves over the data points (not shown). Consequently, one expects that the mean-field approach to calculate the trace of the fabric tensor leads to somewhat overestimated values. For each packing, we calculate the exact value of via Eq. (4) and compare it with the mean-field approximation [Eq. (10)]. Figure 8 reveals that Eq. (10) slightly overestimates in both two- and three-dimensional systems. The deviation increases with the width of the size distribution, but remains less than in all cases. For comparison, note that can reach up to and in 2D and 3D uniform samples, respectively (see Table LABEL:Table2); therefore, ignoring the correction factor would cause up to and error, respectively.

Next, we investigate the average properties of the contact force network. In Sec. III.3, we applied the mean-field approach of Sec. II.3 to estimate the isotropic pressure in a given polydisperse granular sample. However, due to the presence of the normal component of the contact force in Eq. (23), one needs to make one further assumption about the particle-size dependence of to be able to calculate the integral and obtain from the average quantities.

The simulation results [Fig. 9(a)] reveal that the average normal force exerted on the particle is an increasing function of the particle radius for 2D and 3D (the contribution of rattlers is again excluded). With increasing friction coefficient and , the average normal force increases. This is reminiscent of the behavior of as a function of (Figs. 6 and 7). Interestingly, the increasing rates are similar in both figures. Therefore, it is reasonable to assume that the ratio is independent of , as already observed in 2D Madadi05 (). Figure 9(b) confirms the validity of this assumption. We note that the fluctuations in Fig. 9(b) are reduced as the system size increases. In Fig. 10, we compare the exact value of with the corresponding value from Eq. (26) [or Eq. (53)], which is obtained based on the above assumption. The results are in reasonable agreement with theory for both two- and three-dimensional packings, with a standard deviation of to for increasing .

Finally, we turn to the calculation of the stiffness tensor elements for isotropic materials. We note that to evaluate the true elastic moduli, one should apply an incremental strain and measure the resulting change of the stress tensor. Alternatively, one can read the moduli from the elements of the stiffness tensor, assuming the affine motion of the particles, which cannot be taken for granted however, and which is the subject of future studies. Here, using the packing configuration obtained from the simulation, we calculate the elements of the average stiffness tensor via Eq. (31). Next, the elastic moduli of the packing are calculated using Eqs. (32), (33) and (34). The results are then compared to the estimated values of the bulk and shear moduli calculated via Eqs. (42) and (43) [or Eqs. (61) and (62)]. Figure 11 displays the results for several two- and three-dimensional packings; the agreement is satisfactory within a error (also in the case of frictional packings which is not shown here).

According to our analytical results, the ratio between the bulk and shear moduli is for isotropic packings independent of , , and even the size distribution. This suggests that in isotropic packings, the ratio between the -wave velocity and the -wave velocity is always . An experimental test shows that for a compressed polydisperse packing of glass beads remains around over a wide range of pressures from to MPa Lebedev11 () (see also Winkler83 ()). Note, however, that anisotropic regular lattice structures do not necessarily show the same ratio Mouraille06 ().

## Vi Discussion and conclusion

In conclusion, a mean-field approach is developed to isolate the influence of size polydispersity on the physical properties of granular assemblies. We are interested in how the microscale quantities are linked to the macroscale ones.

We find that the trace of fabric and stress tensors factorize into the mean packing properties (for example, average coordination number, packing fraction, and average normal contact force) and dimensionless correction factors, which depend on the moments of the particle-size distribution (and approach unity for monodisperse packings). The method is extended to estimate the elements of the stiffness tensor. This tensor describes the linear affine response of the packing to weak external perturbations, when practically the contact network between the particles remains unchanged. The elements are also proportional to the average quantities and a dimensionless correction factor, which is a function of the size distribution.

Numerical simulations illustrate the validity range of our analytical predictions and of the assumptions on which the mean-field method is based. We note that the deviation of the macroscopic quantities of interest from the average packing properties increases with increasing the width of the particle-size distribution. Figure 12 shows the summarized correction factors as a function of the width of a uniform size distribution, with the average particle size . Neglecting the correction factors would cause remarkable errors, especially for wide distributions. Interestingly, is insensitive to the width of the size distribution in the 3D case. Therefore, according to Eqs. (42) and (43), we expect that the elastic moduli of a polydisperse packing of spheres is only moderately affected by the choice of . The results of molecular dynamics (MD) simulations of soft frictionless spheres imply (see Eq. (12) in Goncu10 ()) that the bulk modulus does not depend on the width of the size distribution, in agreement with our analytical results.

The predictive value of this mean-field method should be examined also by comparing the theoretical predictions with experimental data. For a direct comparison, one needs to measure the average packing properties, e.g.  and , which are not easily accessible in experiments (even though microcomputed tomography (MicroCT) scan determines the geometry with micrometer accuracy nowadays Aste05 ()). Alternatively, by elimination of between our analytical results, one obtains linear relationships between the macroscopic physical properties via some coefficients, which depend on the moments of the size distribution. Such linear relations between macroscopic quantities have been investigated in the literature, e.g. between the elastic moduli and conductivity Bristow60 () or isotropic pressure Mavko03 (), and can be verified experimentally. Future studies will more closely examine the nonaffinity of deformations of isotropic as well as anisotropic packings of frictional and possibly even cohesive particles.

###### Acknowledgements.
We would like to thank T. Unger, L. Brendel, O. Durán, and M. Lebedev for useful discussions and suggestions. M. M. acknowledges financial support by the ANU Digital Core Consortium, and M. R. S. and D. E. W. by DFG Grant No. Wo577/8-1 within the priority program “Particles in Contact”. S. L. acknowledges the support of this project by the Dutch Technology Foundation STW, which is the applied science division of NWO, and by the Stichting voor Fundamenteel Onderzoek der Materie (FOM), financially supported by the Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO).

## Appendix A Analytical results in two dimensions

Fabric tensor . In a two-dimensional packing of disks [Fig. 13(a)], the surface angle covered by a neighboring particle on the reference particle is

 Ω(a)=2arcsin(⟨a⟩a+⟨a⟩), (44)

and the total fraction of the shielded surface is given by

 cs(a)=12πaC(a)∑i=1Ω(a)a=Ω(a)C(a)/2π. (45)

Assuming that is independent of , one can write the mean coordination number as

 z=∫∞0C(a)f(a)da=2πcsq0. (46)

Equations (45) and (46) lead again to Eqs. (9) and (10) for and , with the correction factor

 g1=∫∞0V(a)f(a)Ω(a)daq0∫∞0V(a)f(a)da=⟨a2⟩g⟨a2⟩. (47)

By introducing , we rewrite Eq. (44) as

 Ω(a)=2arcsin(12+ϵ), (48)

and approximate to the first order in for narrow size distributions

 1Ω(a)≃A′1+B′1ϵ, (49)

where and . Figure 13(b) reveals that the approximation has a less than error in the range (or ). Hence, for narrow size distributions becomes

 g1≃1+B′1A′1(⟨a3⟩⟨a⟩⟨a2⟩−1). (50)

Stress tensor . For a two-dimensional disk, by disregarding the direction, i.e., in the plane [by requiring and in Fig. 3(b)], one obtains

 ˜σp=apVp[¯FpnCp∑c=1(cos2(φ)sin(φ)cos(φ)sin(φ)cos(φ)sin2(φ)) (51) +¯Fpt2Cp∑c=1(−sin(φ)cos(φ)cos2(φ)−sin2(φ)sin(φ)cos(φ))],

and its trace

 ˜σpαα=apVpCp∑c=1D∑α=1(¯Fpnnpcαnpcα+¯Fpt2npcαtpc2α) (52)