## Abstract

In this paper, we propose a theoretical model to study the elastic response of a granular material idealized as a random aggregate of identical, elastic, frictional spheres that has been isotropically compressed. When we consider that contacting particles move according to the average deformation, the effective shear modulus is over-predicted with respect to numerical simulation, while the effective bulk modulus is almost captured. We improve upon this simple approach by relaxing the aggregate. The kinematics of a pair of contacting particles is then given by the average deformation and fluctuations in both translations and rotations. We determine analytical expressions for these fluctuations by means of force and moment equilibrium applied to each particle of the pair. In order to derive the incremental stress associated with an incremental deformation, we introduce conditional averages of the fluctuations that are functions of the statistical geometry of the packing. This brings the theoretical predictions of the effective moduli close to those measured in numerical simulations. The variability of the average number of particle contacts per particle is seen to play an important role in the statistical description of the aggregate.

## 1. Introduction

We focus on the initial response of a random aggregate of identical, frictional elastic spheres that has been first isotropically compressed and then incrementally strained. As part of the prediction of the effective moduli, we identify the parameters that characterize the initial state of the aggregate. Although we restrict attention to an initial increment, the corresponding elastic response is important to wave propagation in a granular material that is useful in characterizing the internal structure of the packing (e.g. Domenico 1977; Norris & Johnson 1997; Jia & Mills 2001; Kuwano & Jardine 2002; Agnolin *et al*. 2005).

Digby (1981) and Walton (1987) have obtained expressions for the effective shear modulus and Lame coefficient of such an idealized granular material in terms of the average pressure *p*, the average number of contacts per particle , called the coordination number, the solid volume fraction *v*, and the shear modulus *μ* and Poisson ratio *ν* of the material of the spheres. They assume that the relative displacement of the centres of two contacting particles is given by the average strain and the angular distribution of contacts is isotropic. The effective moduli that they obtain areandwhere(1.1)with the normal component *δ* of the contact displacement given by(1.2)

Jenkins *et al*. (1989) compared these predictions with measurements in numerical simulations of the incremental response of a random aggregate of frictional, elastic spheres. They generated an initial compressed packing at a confining pressure of *p*=138 kPa, with and *ν*=0.63. The measured bulk modulus was slightly different from the average strain prediction; however, the measured shear modulus had a value three times less than that predicted. Another test of the average strain assumption was carried out by Makse *et al*. (1999) who employed an initial state with a similar confining pressure, *p*=100 kPa, and volume fraction, *ν*=0.64, but a relatively high coordination number, . The shear modulus measured in their simulation was approximately 60% lower than that predicted by the average strain theory, while the difference in the measured and predicted bulk modulus was negligible. They also found that the discrepancy between the predicted and measured shear moduli was even more dramatic for frictionless spheres.

Efforts have been made to improve upon the average strain theory. Koenders (1987) introduced fluctuation about the average deformation of contacting particles in a calculation of the effective moduli of idealized assemblies of identical, elastic, frictional discs, but did not make an explicit prediction of effective moduli. Misra & Chang (1993) employed a fluctuation in strain and determined it numerically in the context of a self-consistent scheme for the calculation of the stiffness of the effective media. The calculated stiffness was in good agreement with that measured in numerical simulations on discs, but analytical expressions for the effective moduli were not obtained. Jenkins (1997) and Paine (1997) focused on a typical pair of contacting particles that interact through both the average strain and the fluctuations, assuming that the neighbouring particles moved with the average strain. In this context, Paine (1997) calculated an expression for the effective bulk modulus as a function of a geometrical parameter that she evaluated in a numerical simulation and found that the fluctuations had a small influence. Trentadue (2001) adopted Jenkins's model and, in contrast to what is seen in the simulations, obtained reductions of both the shear and the bulk moduli with respect to the average strain theory by incorporating fluctuations during the initial compression. In addition, Velicky & Caroli (2002) attribute the departure from the pressure dependence of sound waves predicted by the average strain theory to fluctuations of the stress fields. Jenkins *et al*. (2005) applied the pair fluctuation model to an aggregate of frictionless spheres, and found the predicted bulk modulus to be unchanged and the predicted shear modulus reduced to approximately 70% of its average strain prediction. Analytical expressions for the bulk and the shear moduli could be obtained with the pair fluctuation model. Methods that focus on an extended neighbourhood of a single particle, rather than on a typical pair and its neighbours, have also been employed. For example, Jenkins & Koenders (2004), building on the work of Koenders (1994), calculate effective moduli for an isotropic state of frictional, elastic discs that are in reasonable agreement with those measured in numerical simulations by Kuhn (1999).

Here, we extend the calculation of the effective moduli by Jenkins *et al*. (2005) to frictional spheres by introducing fluctuations in both translation and spin to describe the kinematics of a typical pair. The neighbours of the pair are assumed to move with the average strain. We obtain explicit solutions for the fluctuations using force and moment equilibrium. As before, the problem is phrased incrementally, in order that these equations are linear in the fluctuations. The solutions are given in terms of the increment in the average strain and the structural sums that characterize the neighbourhood of the pair. In order to calculate an incremental relation between the average strain and the average stress, we introduce directional averages of these solutions. These involve averages of the structural sums that we express in terms of statistical measures of the neighbourhood of a pair. We find that a measure of particular importance for aggregates with a coordination number less than 6 is the variability of the coordination number. In order to improve the agreement between the predictions and the simulations, we implement the solutions for the fluctuations in the interaction of the pairs with their neighbours and solve the additional fluctuations of the pair. This iteration provides predictions of the effective moduli that are within 10% of those measured in the simulation.

## 2. Theory

### (a) Kinematics

In this section, we attempt to describe the kinematics of an isotropically compressed random aggregate of identical, frictional, elastic spheres. When a subsequent deformation is applied to the aggregate, particles will translate and rotate with respect to each other. In general, the increment in the relative displacement of the contact points of two particles *A* and *B* is given by(2.1)where *d*^{(BA)} is the vector from the centre of *A* to the centre of *B*, and are the increments in translations of the centres of the two spheres, and and are the increments in the rotations about their centres.

Without loss of generality, we focus on average deformations that do not involve average rotations and write the incremental relative displacement of a typical pair of particles in terms of the increment in the average strain and differences and sums of fluctuations(2.2)whereis the increment in the fluctuation in the difference between the translations of the two centres andis the increment in the sum of the fluctuations in the rotations about their centres. The corresponding sum of the fluctuations in the positions of the centres, , and the difference in the fluctuations about the centres, , areandWe note that and do not influence .

Equation (2.1) is a general representation of the kinematics of two contacting particles. However, here we do something simpler and assume that particles in contact with a typical pair translate and rotate with the average motion (Jenkins 1997). If we denote the increment in the translation of centre of the *n*th neighbour of particle *A* by and the increment in its rotation by ,Then, because we assume that for *n*≠*B* only fluctuations in the translation and rotation of particle *A* occur,andTherefore, the incremental contact displacement between particles *n* and *A* is(2.3)

A different approach was proposed by Jenkins & Koenders (2004), who assume that the kinematics of a particle is described by the deformation of the particles in its neighbourhood through a least-squares approximation. In addition, they do not require that the particles of the neighbourhood contact the particle. The kinematics that we employ is somewhat simpler and permits the fluctuations to be determined as solutions of the force and moment balances for the particles of the pair.

### (b) Contact forces

We assume simple relations between the increments in force and contact displacement, and consider each particle to be rigid, with an area of local deformation at each contact.

The incremental contact force exerted by sphere *B* on sphere *A* in a frictional aggregate is not central; there is a normal component in the direction of *d*^{(BA)} and a tangential component that lies in the plane perpendicular to *d*^{(BA)}. This incremental force is usually phrased in terms of the contact stiffnesses and the incremental displacements. Both the normal and the tangential stiffnesses depend upon the compression of the line of centres relative to a configuration in which the distance between the centres is equal to the sum of the radii of the particles. We focus on the first incremental response of the aggregate after an initial isotropic compression. Therefore, it is reasonable to restrict our attention to the elastic response of the medium. Consequently, we write the increment in the contact force in terms of the increment in the relative displacement of the points of contact,(2.4)where ^{(BA)} is the contact stiffness tensor.

The contact stiffness tensor is given in terms of the unit vector by(2.5)where and are the normal and the tangential contact stiffnesses, respectively. The equivalent expression for the incremental contact force between particle *A* and a generic neighbour *n* is(2.6)

## 3. Local equilibrium equations

We have established a simple description of the incremental displacements and contact forces for a typical pair of particles *A* and *B* and their contacting neighbours through the equations (2.2)–(2.4) and (2.6). We next write the equations of force and moment equilibrium for the pair *AB*. We first consider the equations of force equilibrium for particle *A*where *N*^{(A)} is the number of particles in contact with *A*. The corresponding equation for particle *B* isSimilarly, the equations of moment equilibrium for particle *A* arewhile those for particle *B* are

In the equations of equilibrium, we have explicitly distinguished between the contact between the pair and the contacts between a particle of the pair and its other neighbours. However, it is useful to phrase the equilibrium equations in terms of the neighbours of the individual particles of the pair. This will permit the later introduction of an average neighbourhood of a particle. We write, for example,We then introduce tensors that describe the structure of the neighbourhoods of each particle of the pair. In doing this, we distinguish the neighbourhood of particle *A* from that of particle *B* by the order of the indices in the superscriptUsing equation (2.5), we also define three combinations of the components of theseandWe note that is a skew tensor, while is a symmetric tensor and , and, consequently, depend only upon the tangential stiffness.

The equations of force and moment equilibrium for particles *A* and *B* can be written in more compact form using these tensors. For the sake of simplicity, we rewrite only the equations for particle *A*. The equation of its force equilibrium becomes(3.1)while the equation of moment equilibrium is(3.2)

The four equations of equilibrium for the pair *AB*, in matrix form, are(3.3)whereandwithandwhilewith the hat indicating that the first two terms are normalized by the diameter of the particle, and

## 4. Approximate solution

Equation (3.3) represents a system of four equations in the four unknown fluctuations . We note that while the elements in the matrix refer only to the contact *AB*, the components in the matrix take into account the contribution from all of the particles in contact with *A* or *B*. This implies that the tensors in are functions of the number of contacts per particle, *k*. From equation (3.3), we obtainwhere is the identity tensor. An approximate solution can be obtained by neglecting terms proportional to 1/*k*^{2},(4.1)

In order to determine the components of , we evaluate ^{−1} and ^{−1}^{−1}. Here, we note that in equations (3.1) and (3.2), tensors and act on the fluctuations in spin in the balance of forces and the fluctuations in translation in the balance of moment. Since and are odd in ** d**, their averages over the aggregate are 0. Therefore, as an approximation, we may neglect their contribution relative to that of those tensors even in

**. In this case, the problem is uncoupled. In other words, to the order of approximations considered, force equilibrium for a typical particle of a pair depends on the fluctuations in translation, while moment equilibrium is related only to the fluctuations in spin.**

*d*The matrices and becomeandAs can easily be verified,andwhereWith these, the solutions for the fluctuations are(4.2)(4.3)and(4.4)(4.5)The kinematics of the pair fluctuation model is completed by these terms. We next turn to the determination of the average stress and its evaluation in terms of averages of the local fluctuations.

## 5. Stress tensor

The expression for the contact force between a pair of particles can be employed to evaluate the average stress tensor of the aggregate. The general expression for the incremental average stress tensor for an aggregate of *N* particles is (e.g. Love 1927, note B)(5.1)where *V*^{p} is the volume of the cell that contains particle *p* in a Voroni tessellation of the aggregate. For simplicity, we associate with each particle the average cell volume *V*,Then, equation (5.1) can be written as(5.2)where *q* is the total number of contacts in the aggregate.

If we had detailed knowledge of the tensors , , and for each pair in the aggregate, we could evaluate equations (4.2) and (4.5) for each pair and, using equations (2.2) and (2.4), calculate the average stress tensor by means of equation (5.2). However, such detailed information is available only from the numerical simulations. To overcome this problem, we first rewrite equation (5.2) by collecting the contact forces for pairs of particles whose contact vectors are close to a given direction. We take *M*^{p} to be the number of pairs with contact vectors within the *p*th element of *R* equal elements Δ*Ω*^{p} of solid angle in which the surface of the unit sphere has been partitioned such thatThen

Next, following Jenkins *et al*. (2005), we introduce a conditional average aswhere *F*^{(CD)} is a contact force for a pair within the element of solid angle centred on . Then, upon assuming that ** F** and

**are not correlated and with , the incremental stress tensor becomes(5.3)The corresponding conditional average of the incremental contact force isand equation (5.3) becomeswhere we recall that**

*d**N*is the total number of particles in the aggregate, equal to twice the number of pairs and

*V*

^{tot}is the total volume of the aggregate.

In the limit that Δ*Ω* tends to d*Ω*(5.4)where *n*=6*v*/*πd*^{3} is the number of particles per unit volume and is the number of contacts in the element d*Ω* of solid angle centred at . We also have thatEquation (5.4) is the continuous form for the increment in the average stress tensor. Its value depends upon the determination of the conditional average . Here, we replace the generic direction with and, using equation (2.4), write(5.5)

Both the stiffness and the displacement can be expressed as the sum of a conditional average and a fluctuation. We assume that the contact stiffness, the contact vector of the pair and the fluctuations are not correlated within a given element of solid angle. Consequently,(5.6)Upon taking the conditional average of equations (4.2) and (4.5), we obtain(5.7)and(5.8)Since the existence of the contact between particles *A* and *B* provides a symmetry, on an average, about the plane perpendicular to , and . This is because the representation for is function of an odd number of unit vectors ; while, for , the number of unit vectors is even. We also have, for example, and . Furthermore, because contributions to the fluctuations related to ^{−1}^{−1} are 1/*k* lower than those related to ^{−1}, we adopt for the former simple conditional averages and, for the inverses, assume that the conditional average of the inverse is the inverse of the conditional average. Consequently, upon taking the conditional averages of equations (4.3) and (4.4),At this point, we have all the elements to carry out the integral in equation (5.4).

## 6. Averages

In the expression (5.6), we need to evaluate , and . To do so, we first assume that the conditional averages of the stiffnesses are equal to their averages over all the contacts of the aggregate: and . Our justification for this is that the packing is isotropic and subjected to an isotropic compression. We identify the average values and with those determined using the average strain assumption. In other words, they are determined through equation (1.1), where we adopt the average strain relation proposed by Jenkins *et al*. (1989) between the normal component of the contact displacement *δ* and the average pressure *p*, equation (1.2).

Since we deal with an aggregate of identical spheres for which the ratio between normal and tangential contact stiffnesses is close to unity, we neglect terms in the fluctuation that are proportional to the difference in the contact stiffness. In this case,is adopted andThen, from equation (5.7),(6.1)From equation (5.8) and the definitions of the tensors , and The tensors and have the same symmetries and can be represented by a fourth rank tensor ,(6.2)Then,(6.3)where only the coefficients *ξ*_{3} and *ξ*_{5} survive the contraction with the alternating tensor. Since we are interested in the first incremental response of the aggregate, after the material has been isotropically compressed, the coefficients *ρ* and *ξ* in equations (6.1) and (6.2) are functions only of the statistics of the particles in contact with the pair *AB*. In an initial state that involved deviatoric strain, these coefficients would also be a function of the applied strain.

## 7. Elastic moduli

With the incremental contact force given by equation (5.6), and the expressions for the fluctuations (6.1) and (6.3), we employ the identitiesandand calculate the increment in average stressConsequently, the effective elastic moduli areandThe bulk modulus, calculated in terms of these, isObviously, the bulk modulus is not influenced by the tangential stiffness or, equivalently, by the fluctuations in spin.

## 8. Statistical models

We focus our attention on equations (6.1) and (6.3) and attempt to model these averages. Within a given element of solid angle, each tensor can be represented as the sum of average and fluctuation. The fluctuations occur because the number and direction of contacts will change from particle to particle for the pairs in the element. However, owing to the symmetry of the local packing, those tensors whose representation involves an even number of vectors *d*^{(BA)} will fluctuate less than those that involve an odd number. For tensors like , and , we take this into account in the average of the product of fluctuations by neglecting the variability in the direction of the contacts and retaining only the variability associated with the number of contacts.

We denote the fluctuation of tensor by ′. The inverse of can be written in terms of its average and fluctuation byup to an error proportional to the cube of the fluctuation. Then, in the equation (5.7), the first contribution becomes(8.1)Similarly, for ** B**,so in equation (5.8), the first contribution becomes(8.2)

### (a) Simple conditional average

In replacing the summation with an integral, we incorporate distributions of both the number of contacts per particle and the orientation of the contacts. Then, for particles for which the probable number of contacts is *f*(*k*)d*k*, *g*(α)d*Ω* is the number of these within the solid angle d*Ω* centred at *α*. We assume that *α* and *k* are two random independent variables with the joint probability *l*(α,*k*)=*g*(α)*f*(*k*), defined so that *l*(α,*k*)d*Ω* d*k* is the number of contacts in an element of solid angle d*Ω*=sin *θ* d*θ* d*ϕ* centred at *α*, given that there is a contact at .

The presence of particle *B*, located at the pole in the upper hemisphere excludes other particles from the interval 0≤*θ*≤*π*/3 and influences the distribution of particles at orientations for *θ*>*π*/3 in a way that decreases with distance from the pole. In order to model this in a rough way, we assume that the contacts are uniformly distributed within three intervals, so that within *π*/3≤*θ*≤*π*/2, there are (2/3)(*k*/2−1+*k*/4) contacts; within *π*/2≤*θ*≤2*π*/3, there are (1/3)(*k*/2−1+*k*/4) contacts and within 2*π*/3≤*θ*≤*π*, there are *k*/4 contacts. This implies that the last interval is independent of the presence of particle *B* in upper hemisphere. This is a slightly more accurate distribution than that employed by Jenkins *et al*. (2005). With these hypotheses, the distribution function is

Then, for example, the continuous version for the conditional average of *A*^{(BA)} iswhere the integration is over all solid angles *Ω*_{α} consistent with the presence of particle *B*, and over all possible numbers of contacts per particles. The first term is associated with the presence of the particle *B*. When expressed in terms of the distribution functions *f*(*k*) and *g*(*α*), this average iswhere by definition,

It is a straightforward calculation to determine thatandwhere , in which the overbar indicates that the coefficient depends upon . We also havewhere and . Then, given and ,(8.3)where . Its inverse can be approximated by taking for its isotropic partwith . Finally,

### (b) Average of the product of the fluctuations

We next evaluate the second terms on the right side of equations (8.1) and (8.2). By definition,wherewithWe note that we now must consider an average of a double summation related to the same particle. The continuous form of this is (Jenkins *et al*. 2005)By virtue of the expressions of *g*(*α*),We write the coordination number as the sum of its average value and a fluctuationThen, upon evaluating the integral,In addition,andFinally, equation (8.1) may be written as

In equation (8.2), by definition,whereSince and are tensors of even rank, we neglect fluctuations in orientation and retain fluctuation in the number of contacts per particle. To do the latter, we extend the formulation of Jenkins *et al*. (2005) to include *k*, a random variable independent for the orientations *α* and *β*, and write(8.4)where *F*(*α*, *β*) is the joint probability density, defined so that *F*(*α*, *β*)d*Ω*_{α}d*Ω*_{β} is the fraction of contacts with *α* in d*Ω*_{α} and *β* in d*Ω*_{β}, including the possibility that the two directions can coincide.

At this point, we simplify the problem by assuming that *α* and *β* are independent and write equation (8.4) asThe final result can be represented bywhereandFocusing on ,(8.5)

## 9. Prediction

With the evaluation of the coefficients of the tensors introduced in the expression for the fluctuations in spin and translation, we now have all of the elements to predict the effective elastic moduli. In particular, with respect to equation (6.1),(9.1)while for the tensor , according to equations (6.2) and (8.2),(9.2)and *ξ*_{3}=0. Both *ρ* and *ξ*_{5} contain three contributions. The first two are related to the expression for ^{−1} given in equation (4.1). In particular, the second contribution incorporates fluctuations in the structure that, in this case, are only functions of the fluctuations in the coordination number. The last contributions to *ρ* and *ξ*_{5} are terms proportional to 1/*k* in equation (4.1) and, therefore, their contribution increases for decreasing values of *k*. We also note that in equation (8.5), the terms within the parenthesis on the right-hand side are the elements of the representation for a tensor like in equation (6.2). Finally, the hypothesis that depends only upon the fluctuation in the number of contacts per particle, leads to a simpler representation with *ξ*_{3}=0.

Therefore, the solutions for the fluctuations are(9.3)and(9.4)Consequently, the corresponding expressions for effective elastic moduli arewhere , andWith these, the effective bulk modulus isThese theoretical expressions incorporate contributions from the average and fluctuations. We note that the dependence on the ratio of the mean square of *k*′ and the square of . Such a term also appears in Jenkins & Koenders (2004), where the authors emphasize its importance to the prediction of the stress–strain behaviour of an aggregate of frictional discs with an average coordination number of approximately 4.

## 10. Numerical simulation

In order to test the theoretical prediction, we consider a numerical simulations carried out using the code TRUBAL developed by Cundall (1988). Each simulation employed spheres of equal radii: *R*=0.1×10^{−3} m. The shear modulus of the material of the spheres was *μ*=2.9×10^{10} Pa and the Poisson ratio was *ν*=0.2. This simulation employed a system of 10 000 spheres. The initial state was obtained in the manner described by Makse *et al*. (1999). In particular, we take two initial isotropic states characterized by the same isotropic pressure, *p*=100×10^{3} Pa, a slight variation of the solid volume fraction approximately *v*=0.638, and with two different coordination numbers, and and, respectively, with in the former case and in the latter. For these packings, the number of rattlers (particles with no force) is small, so only slight variations in and *k*′ result from ignoring them; however, the influence of the rattlers on the coordination number is more pronounced for when . We also note that the ratio of, for example, the normal contact stiffness measured in the simulation and that given by equation (1.1) varies from 0.97 when to 0.92 when These states, for random packings of spheres, are close to what obtained by Agnolin *et al*. (2005) where they create similar dense aggregates, under the same pressure, with different coordination numbers.

The results for the elastic moduli, in mega Pascal, are given by table 1 and table 2. This shows the decrease of the predicted effective moduli from the average strain assumption. Moreover, it seems that the variability in the number of contacts per particle plays a more important role for the second packing rather that in the case of . However, we note that a discrepancy between theory and simulation is still rather large. We attempt to reduce this gap by incorporating more relaxation in the kinematics by introducing an iteration method based upon the pair fluctuation solution.

## 11. Iteration

In order to improve on the pair fluctuations, we introduce a model in which the kinematics of a particle and its contacting neighbours is given by (2.2). We assume that in a typical neighbour in the aggregate each pair moves according to the average deformation and a fluctuation and that each pair is independent of the others. We then focus on a typical pair with a definite orientation. In order to guarantee the equilibrium of this pair, we introduce an additional fluctuation and assume that its neighbours move according to equation (2.2), i.e. for a pair *AB*,where is given by the conditional average of (2.2) and the additional fluctuations are labelled with the index 1. For the *n*th particle in contact with *A*,

The equations of force and moment equilibrium for particle *A*, when phrased in terms of the new unknowns, have the same structure of equations (3.1) and (3.2) if, for example, we replace by , ^{(BA)} bywhereand ^{(BA)} byThe second term on the right-hand side of the expression for ^{(1)}*J*^{(BA)} may be neglected with respect to the first, as both *ξ*_{5} and *ρ* are much smaller than unity.

The solutions are easily derived by neglecting terms, in this case, proportional to 1/*k*and

If we again take the conditional averages,andwith . The additional fluctuations have the representationswhere , andwhere .

With these, the effective shear and bulk moduli areandwhere again

The results for the effective elastic moduli, in mega Pascal, can be collected in table 3 and table 4.

## 12. Conclusion

We have derived expressions for the effective elastic moduli for an isotropic random aggregate of identical, elastic, frictional spheres. The theory is based upon a plausible kinematics for the contacting particles described by averages and deviations from them. We have obtained analytical expressions for both the fluctuations in translation and rotation using force and moment equilibrium applied to each particle of a pair. In order to evaluate the influence of the fluctuations on the elastic response of the aggregate, we have introduced conditional averages for them that depend on the geometry of the packing. In particular, the variation of the number of contacts per particle can play an important role. Here, we have attempted to derive simple expressions for the effective moduli of an isotropic, random aggregate of identical, frictional, elastic spheres. The derivation has involved numerous approximations, some of them uncontrolled. These include the initial assumption of pair fluctuations, the neglect of terms proportional in the solution of the force and moment equilibrium equations, the neglect of the interaction of the spin and translational fluctuations in the equilibrium equations, and the neglect of terms proportional to the cube of the fluctuations in the structural tensors. In making these approximations, we have been guided by numerical simulations that will be reported on in detail elsewhere. We believe that the elements that we have included in our derivation are those that are essential to accurate approximate prediction of the elastic moduli. The predicted moduli are in excellent agreement with those measured in numerical simulations.

## Acknowledgments

This research was supported by M.I.U.R.-Cofin2005, Ministero Affari Esteri, Direzione Generale per la Promozione e la Cooperazione Culturale. The authors are grateful to Vanessa Magnanimo for providing us with the results of her numerical simulations.

## Footnotes

- Received August 28, 2006.
- Accepted November 10, 2006.

- © 2006 The Royal Society