# Why Effective Medium Theory Fails in Granular Materials

## Abstract

Experimentally it is known that the bulk modulus, $`K`$, and shear modulus, $`\mu `$, of a granular assembly of elastic spheres increase with pressure, $`p`$, faster than the $`p^{1/3}`$ law predicted by effective medium theory (EMT) based on Hertz-Mindlin contact forces. Further, the ratio $`K/\mu `$ is found to be roughly pressure independent but the measured values are considerably larger than the EMT predictions. To understand the origin of these discrepancies, we have undertaken numerical simulations of a granular assembly of spherical elastic grains. Our results for $`K(p)`$ and $`\mu (p)`$ are in good agreement with the existing experimental data. We show, also, that EMT can describe their pressure dependence if one takes into account the fact that the number of grain-grain contacts increases with $`p`$. Most important, the affine assumption (which underlies EMT), is found to be valid for $`K(p)`$ but to breakdown seriously for $`\mu (p)`$. This explains why the experimental and numerical values of $`\mu (p)`$ are much smaller than the EMT predictions.

PACS: 81.06.Rm, 81.40.Jj

The study of sound propagation and nonlinear elasticity in unconsolidated granular matter is a topic of great current interest . In the simplest experiments, a packing of glass beads is confined under hydrostatic conditions and the compressional and shear sound speeds, $`v_p`$ and $`v_s`$, are measured as functions of pressure, $`p`$ . In the long-wavelength limit, the sound speeds are related to the elastic constants of the aggregate: $`v_p=\sqrt{(K+4/3\mu )/\rho ^{}}`$ and $`v_s=\sqrt{\mu /\rho ^{}}`$, where $`\rho ^{}`$ is the system’s density. In a recent Letter  acoustic measurements were made on bead packs under uniaxial stress and it was suggested that long wavelength compresional waves can be described in terms of an effective medium. Thus, it would of great value to have a reliable EMT to describe sound propagation as a function of applied stress. However, our analysis, together with the work of others, raises serious question about the validity of the generally accepted theoretical formulation. The EMT  predicts that $`K`$ and $`\mu `$ both vary as $`p^{1/3}`$, and that the ratio $`K/\mu `$ is a constant (independent of pressure and coordination number) dependent only on the Poisson’s ratio of the material from which the individual grains are made.

Experimentally (see Fig. 1), the bulk and shear moduli increase more rapidly than $`p^{1/3}`$ and the values of $`K/\mu `$ are considerably larger than the EMT prediction. These discrepancies between theory and experiment could be due to the breakdown of the Hertz-Mindlin force law at each grain contact as proposed in  for the case of metallic beads with an oxide layer, and in  for grains with sharp angularities. Alternatively, they could be associated with the breakdown of some of the assumptions underlying the EMT, for example, that the number of contacts per grain is pressure independent, which may not be the case as several authors have suggested .

In this Letter we report calculations of $`K(p)`$ and $`\mu (p)`$ based on granular dynamics (GD) simulations using the Discrete Element Method developed by Cundall and Strack . Here it is assumed at the outset that one has an assembly of spherical soft grains which interact via the Hertz-Mindlin force laws. We find good agreement with the existing experimental data, thus confirming the validity of Hertz-Mindlin contact theory to glass bead aggregates. Further, we can explain the two problems with EMT described above. First, if the calculated increase of the average coordination number with $`p`$ is taken into account, the modified EMT gives an accurate description of the bulk modulus found in the simulations, $`K(p)`$; for $`\mu (p)`$ we obtain a curve whose shape is in good agreement with the simulation data but whose values are seriously offset therefrom. Second, the EMT makes the affine assumption in which the motion of each grain is specified simply in terms of the externally applied strain (see below). We show that while the affine assumption is approximately valid for the bulk modulus, it is seriously in error for the shear modulus; this is why the EMT prediction of $`K/\mu `$ differs significantly from the experimental value.

Numerical Simulations. At the microscopic level the grains interact with one another via (1) non-linear Hertz normal forces and (2) friction generated transverse forces. The normal force, $`f_n`$, has the typical 3/2 power law dependence on the overlap between two spheres in contact, while the transverse force, $`f_t`$ depends on both the shear and normal displacements between the spheres . For two spherical grains with radii $`R_1`$ and $`R_2`$:

$$f_n=\frac{2}{3}C_nR^{1/2}w^{3/2},$$
(2)   
$$\mathrm{\Delta }f_t=C_t(Rw)^{1/2}\mathrm{\Delta }s.$$
(3)   
Here $`R=2R_1R_2/(R_1+R_2)`$, the normal overlap is $`w=(1/2)[(R_1+R_2)|\stackrel{}{x}_1\stackrel{}{x}_2|]>0`$, and $`\stackrel{}{x}_1`$, $`\stackrel{}{x}_2`$ are the positions of the grain centers. The normal force acts only in compression, $`f_n=0`$ when $`w<0`$. The variable $`s`$ is defined such that the relative shear displacement between the two grain centers is $`2s`$. The prefactors $`C_n=4G/(1\nu )`$ and $`C_t=8G/(2\nu )`$ are defined in terms of the shear modulus $`G`$ and the Poisson’s ratio $`\nu `$ of the material from which the grains are made. In our simulations we set $`G=29`$ GPa and $`\nu =0.2`$. We assume a distribution of grain radii in which $`R_1=0.105`$ mm for half the grains and $`R_2=0.095`$ mm for the other half. Our results are, in fact, quite insensitive to the choice of distribution. We also include a viscous damping term to allow the system to relax toward static equilibrium .

Our calculations begin with a numerical protocol designed to mimic the experimental procedure used to prepare dense packed granular materials. In the experiments the initial bead pack is subjected to mechanical tapping and ultrasonic vibration in order to increase the solid phase volume fraction, $`\varphi _s`$. The simulations begin with a gas of 10000 spherical particles located at random positions in a periodically repeated cubic unit cell approximately 4 mm on a side. At the outset, the transverse force between the grains is turned off ($`C_t=0`$). The system is them compressed slowly until a specified value of $`\varphi _s`$ is attained (see dashed lines in Fig. 2). The compression is then stopped and the grains are allowed to relax. If the compression is stopped before reaching the critical volume fraction, $`\varphi _s0.64`$, corresponding to random close packing (RCP), the system will relax to zero pressure and zero coordination number, since the system cannot equilibrate below RCP. The compression is then continued to a point above the critical packing fraction and the target pressure is maintained with a “servo” mechanism  which constantly adjusts the applied strain until the system reaches equilibrium. Because there are no transverse forces, the grains slide without resistance during the relaxation process and the system reaches the high volume fractions found experimentally.

The simulated granular aggregate relaxes to equilibrium states in which the average coordination number, $`<Z(p)>`$, increases with pressure as seen in Fig. 2. For low pressures compared with the shear modulus of the beads $`<Z>6`$, while in two dimensions the same preparation protocol gives $`<Z>4`$. Such low coordination numbers can be understood in terms of a constraint argument for frictionless rigid balls , which gives $`<Z>=2D`$, where D is the dimension. These values should be valid in the limit of low pressure when the beads are minimally connected near RCP  (or in the rigid ball limit $`G\mathrm{}`$). For large values of the confining pressure more grains are brought into contact, and the coordination number increases . Empirically, we find that

$$<Z(p)>=6+\left(\frac{p}{0.06\text{MPa}}\right)^{1/3}.$$
(4)   
Comparison with Experiment. Consider, now, the calculation of the elastic moduli of the system as function of pressure. Beginning with the equilibrium state describe above, we first restore the transverse component of the contact force interaction. We then apply a small perturbation to the system and measure the resulting response. The shear modulus is calculated in two ways, from a pure shear test, $`\mu =(1/2)\mathrm{\Delta }\sigma _{12}/\mathrm{\Delta }ϵ_{12}`$, and also from a biaxial test, $`\mu =(\mathrm{\Delta }\sigma _{22}\mathrm{\Delta }\sigma _{11})/2(\mathrm{\Delta }ϵ_{22}\mathrm{\Delta }ϵ_{11})`$. The bulk modulus is obtain from a uniaxial compression test, $`K+4/3\mu =\mathrm{\Delta }\sigma _{11}/\mathrm{\Delta }ϵ_{11}`$. Here the stress, $`\sigma _{ij}`$, is determined from the measured forces on the grains , and the strain, $`ϵ_{ij}`$, is determined from the imposed dimensions of the unit cell.

In Fig. 1 our calculated values of the elastic moduli as a function of pressure are compared with EMT and with experimental data. Because there is a considerable degree of scatter in the experimental results we performed our own experiments with standard sound propagation techniques. A set of high quality glass beads of diameter $`45\mu `$m are deposited in a flexible container of 3 cm height and 2.5 cm radius. Transducers and a pair of linear variable differential transformers (for measurement of displacement) were placed at the top and bottom of the flexible membrane, and the entire system was put into a pressure vessel filled with oil. Before starting the measurements, a series of tapping and ultrasonic vibrations were applied to the container in order to let the grains settle into the best possible packing. We then applied confining pressures ranging from 5 MPa to 100 MPa. The pressure was cycled up and down several times until the system exhibited minimal hysteresis. At this point shear and compressional waves were propagated by applying pulses with center frequencies of 500 KHz. The sound speeds and corresponding moduli were obtained by measuring the arrival time for the two sound waves.

From Fig. 1 we see that our experimental and numerical results are in reasonably good agreement. Also shown are measured data from Domenico  and Yin . Clearly, the experimental data are somewhat scattered. This scatter reflects the difficulty of the measurements, especially at the lowest pressures where there is significant signal loss. Nevertheless, our calculated results pass through the collection of available data. Also shown in Fig. 1 are the EMT predictions 

$$K=\frac{C_n}{12\pi }(\varphi _sZ)^{2/3}\left(\frac{6\pi p}{C_n}\right)^{1/3},$$
(6)   
$$\mu =\frac{C_n+(3/2)C_t}{20\pi }(\varphi _sZ)^{2/3}\left(\frac{6\pi p}{C_n}\right)^{1/3}.$$
(7)   
The EMT curves are obtained using the same parameters as in the simulations; we also set $`Z=6`$ and $`\varphi _s=0.64`$, independent of pressure. At low pressures we see that $`K`$ is well described by EMT. At larger pressures, however, the experimental and numerical values of $`K`$ grow faster than the $`p^{1/3}`$ law predicted by EMT. The situation with the shear modulus is even less satisfactory. EMT overestimates $`\mu (p)`$ at low pressures but, again, underestimates the increase in $`\mu (p)`$ with pressure.

To investigate the failure of EMT in predicting the correct pressure dependence of the moduli, we plot in Fig. 3 the shear modulus divided by $`p^{1/3}`$. For such a plot, EMT predicts a horizontal straight line but we see that the numerical and experimental results are clearly increasing with $`p`$. We can explain this behavior by modifying Eq. (3b) to take account of the pressure dependence of the coordination number $`<Z(p)>`$ from Fig. 2 and the pressure dependence of $`\varphi _s(p)`$ (which is a much smaller effect). The modified EMT is also plotted in Fig. 3 and we see that it predicts the same trend with pressure as the simulations. The experimental data also seem to be following this trend but more data over a larger pressure range are clearly needed. Not shown in Fig. 3 is a similar analysis of $`K(p)`$ but the result is that the modified EMT is in essentially exact agreement with our numerical simulations. It is for this reason that we focus on $`\mu (p)`$.

Another way of seeing the breakdown of EMT is to focus on the ratio $`K/\mu `$. According to Eqs. (7), $`K/\mu =5/3(2\nu )/(54\nu )`$, independent of pressure, a value which depends only on the Poisson’s ratio of the bead material. The experiments give $`K/\mu 1.11.3`$. Our simulations give $`K/\mu 1.05\pm 0.1`$ in good agreement with experiments. EMT predicts $`K/\mu =0.71`$, if we take $`\nu =0.2`$ for the Poisson’s ratio of glass. \[The EMT prediction is quite insensitive to variations of $`\nu `$; $`K/\mu =0.71\pm 0.04`$ for $`\nu =0.2\pm 0.1`$.\]

To understand why $`\mu `$ is over predicted by EMT we must examine the role of transverse forces and rotations in the relaxation of the grains. \[These effects do not play any role in the calculation of the bulk modulus.\] Suppose we re-define the transverse force by introducing a multiplicative coefficient $`\alpha `$, viz: $`\mathrm{\Delta }f_t=\alpha C_t(Rw)^{1/2}\mathrm{\Delta }s`$; with $`\alpha =1`$ we recover our previous results. To quantify the role of the transverse force on the elastic moduli, we calculate $`K(\alpha )`$ and $`\mu (\alpha )`$ at a given pressure $`p=100`$ KPa \[Fig. 4a\]. This pressure is low enough that the changing number of contacts is not an issue. Surprisingly, the shear modulus becomes negligibly small as $`\alpha 0`$. As expected, $`K`$ is independent of the strength of the transverse force. To compare with the theory we also plot the prediction of the EMT, Eq. (7) in which $`C_t`$ is rescaled by $`\alpha C_t`$. We see that the EMT fails in taking into account the vanishing of the shear modulus as $`\alpha 0`$. However it accurately predicts the value of the bulk modulus, which is independent of $`\alpha `$.

There are two main approximations in the EMT: (1) All the spheres are statistically the same, and it is assumed that there is an isotropic distribution of contacts around a given sphere. (2) An affine approximation is used, i.e., the spheres at position $`X_j`$ are moved a distance $`\delta u_i`$ in a time interval $`\delta t`$ according to the macroscopic strain rate $`\dot{ϵ}_{ij}`$ by $`\delta u_i=\dot{ϵ}_{ij}X_j\delta t`$. The grains are always at equilibrium due to the assumption of isotropic distribution of contacts and further relaxation is not required.

In the GD calculation of the shear modulus an affine perturbation is first applied to the system. The shear stress increases (from A to B in Fig. 4b) and the grains are far from equilibrium since the system is disordered. The grains then relax towards equilibrium (from B to C), and we measure the resulting change in stress from which the modulus is calculated. To better understand the approximations involved in the EMT, suppose we repeat the GD calculations taking into account only the affine motion of the grains and ignoring the subsequent relaxation. The resulting values of the moduli are plotted in Fig. 4a as open symbols and we see that the moduli calculated this way are very close to the EMT predictions. Thus, the difference between the GD and EMT results for the shear modulus lies in the non-affine relaxation of the grains; this difference being largest when there is no transverse force. By contrast, grain relaxation after an applied isotropic affine perturbation is not particularly significant and the EMT predictions for the bulk modulus are quite accurate.

The surprisingly small values of $`\mu `$ found as $`\alpha 0`$ can be understood as a melting of the system that occurs when the system is close to the RCP fraction. This fluid like behavior (when $`C_t=0`$) is closely related to the melting transition seen in compressed emulsions  and foams . At the RCP fraction the system behaves like a fluid with no resistance to shear. By contrast, molecular dynamics simulations of glasses, in which the atoms interact by purely longitudinal forces, predict non-vanishing shear speeds . The crucial difference between these two systems is the local coordination of the particles. In the granular system, the coordination number near RCP (where the balls are only weakly deformed) is $`<Z>=6`$; the system is at its minimal coordination number. In glasses, however, the number of neighbors is closer to 10 and the motion of the grain is highly constrained.

Conclusions. Our GD simulations are in good agreement with the available experimental data on the pressure dependence of the elastic moduli of granular packings. They also serve to clarify the deficiencies of EMT. Grain relaxation after an infinitesimal affine strain transformation is an essential component of the shear (but not the bulk) modulus. This relaxation is not taken into account in the EMT. In the limit $`\alpha 0`$ a packing of nearly rigid particles responds to an external isotropic load with an elastic deformation and a finite $`K`$. By contrast, such a system cannot support a shear load ($`\mu 0`$) without severe particle rearrangements. This may indicate a “fragile” state of the system  where inter particle forces are organized along “force chains” (stress paths carrying most of the forces in the system) oriented along the principal stress axes. Such fragile networks support, elastically, only perturbations compatible with the structure of force chains and deform plastically otherwise. Clearly, there is a need for an improved EMT; recent work on stress transmission in minimally connected networks  may provide an alternative formulation and allow to describe properly the response of granular materials to perturbations.

ACKNOWLEDGMENTS. We would like to thank J. Dvorkin, J. St. Germain, B. Halperin, J. Jenkins, and D. Pissarenko for many stimulating discussions.