## Abstract

We consider a random aggregate of identical, frictionless spheres whose contact is maintained by an applied pressure. The aggregate is then subjected to an axial compression at fixed pressure. We show that the incremental elastic response of the resulting transversely isotropic material is characterized by six rather than by five independent coefficients and that the stiffness tensor does not have the major symmetry. This is because we permit deviations from an affine deformation that are determined by local equilibrium, when anisotropy is present. Discrete element numerical simulations confirm these findings.

## 1. Introduction

In his treatise on the theory of elasticity [1], Love dedicates a section in the appendix (note B) to address some issues on the Cauchy molecular theory (also called average strain or affine deformation theory) and its relation with the continuum approach based on the existence of a potential. The Cauchy model predicts, in the isotropic case, one independent modulus for the aggregate, whereas the continuum model has two independent moduli, as expected. The interesting point is that the appendix is not only a simple historical review of a famous dispute between prominent scientists of the nineteenth century (see also [2]), but it provides valuable comments. In particular, Love underlines that the simple symmetric arrangement of the ‘structure elements’ and the absence of equilibrium equations explicitly treated are two key points in understanding the weakness of the molecular model. This model, in addition to its historical role, has been adopted more recently to predict the response of aggregate of particles like granular materials [3]. The result is interesting because it permits to understand the macroresponse of the aggregate from the local particles interaction. On the other hand, the comparison of the affine deformation theory against numerical simulation and physical experiments has shown a weakness in the model prediction. Here, we focus on a simpler aggregate than a granular material, made of frictionless, random, elastic particles [4]. We do this in the context of a theory that goes beyond the simple average strain. An example of such an aggregate are compressed emulsions, at rather high volume fraction, made by immiscible mixture of two fluids, one dispersed in the other [5,6]. When a surfactant is introduced to stabilize the interfaces of the droplets, the emulsions are seen as a collection of deformable, frictionless particles interacting only through normal forces [7]. Unlike granular materials, particles in emulsions do not transmit a tangential force but, like granular materials, emulsions have a clear and well-defined elastic response [5,8].

We first show the prediction of the elastic response of an isotropic and transversely isotropic aggregate deforming with an affine strain. We have, respectively, one and three independent moduli, instead of the expected two and five. Next, we extend the model and, following Jenkins *et al.* [9], assume that contacting particles move with average deformation and a fluctuation. We explicitly employ equilibrium equations, unlike the Cauchy model, and we determine the fluctuations. We next introduce suitable averages for the fluctuations to calculate an incremental stress–strain relation. In the isotropic case, it is known that the fluctuation model predicts two independent moduli [9], as seen in simulation and physical experiments. Seen from the historical point of view, this result suggests that it is not necessary to introduce a non-central force (in Voigt [10] following Poisson [11] in their study on crystal models) to derive two independent moduli.

In the case of a transversely isotropic aggregate, the present theory predicts six independent moduli. In particular, the fluctuation model predicts an elastic tensor without major symmetry and, therefore, it is not possible to define a potential for the aggregate. This is the relevant result of the paper that is obtained, because particles can deviate from their average motion in the presence of anisotropy. That is, while previously fluctuations have been introduced to obtain a better quantitative prediction than in the simple average theory [9,12–14], here, we show that stiffness tensor **A** of a transversely isotropic, assembly of frictionless, elastic particles is characterized by the lack of major symmetry, say *A*_{ijkl}≠*A*_{klij}. Numerical simulations are carried out to test the theory and indicate the absence of the major symmetry of the stiffness tensor in a transversely isotropic aggregate of frictionless particles.

## 2. Average strain

The incremental response of a random, frictionless aggregate of *N* identical, elastic spheres with diameter *D*, first isotropically compressed and then sheared, is derived. We first review the molecular Cauchy approach, which is equivalent to assuming that the incremental, relative displacement of the centres of contacting particles, *K*_{N} is the normal contact stiffness

Here, *G* and *ν* are, respectively, the shear modulus and the Poisson ratio of the material of the particle, whereas *δ* is the compressive displacement of the centres of the particles:
*E*_{ql}=−(*ε*/3)*δ*_{ql}, where *δ*_{ql} is the Kronecker's delta, *K*_{N} is constant and independent of *ε*=−*E*_{kk} (positive in a decrease in volume) in the initial compressed state. The analytical expression for the stress increment ([1], note B) is written in terms of the number of particles per unit volume *n* and a contact distribution function *k* is the coordination number (the average number of contacts per particle). Then
**A,** *A*_{ijkl}=*A*_{klij}, is obtained.

We now consider a uniaxial compression, with **h**≡**y**_{3} be the axis of compression; the average strain is
*ε* is the total volume strain including the part associated with the initial isotropic compression. We restrict our analysis to the range of deformation in which the deviatoric part of the strain is small compared with the isotropic part. This is a rather small regime of deformation, for which anisotropy induced by strain is already present. In this range, it is also plausible to assume no change in the geometric contact network, the contact distribution is still isotropic. Both assumptions, anisotropy and small variation of fabric, are tested in the simulation, whose results we show later in the paper. A detailed numerical analysis of fabric and anisotropy can be found, for example, in references [15–17], in which the range of deformation of interest is quite large.

Given equations (2.3), (2.4) and (2.9), we note that the stiffness varies with the orientation of the contact with respect to **h**, and this is the way induced anisotropy enters in the problem. We expand the term in equation (2.3) associated with *δ* and retain the contributions that are linear in the ratio of the deviatoric to the isotropic strain:
*ρ*=*nD*^{3}*G*/[2(1−*ν*)3^{1/2}]. The general representation for the elastic tensor is
*et al.* [18], where anisotropy is associated with the particle's arrangement. Here, we differ as the aggregate is frictionless, and the anisotropy is induced rather than inherent.

We attempt to improve on the average strain theory by introducing fluctuations, following a theoretical approach adopted to predict the behaviour of a granular material [9,12,19]. In §3, we elaborate upon this more sophisticated theory; equilibrium is treated explicitly and a possible link between the deformation of the pair and its neighbours is introduced. The result is that, in the isotropic case, two independent effective moduli, in agreement with both numerical simulation and physical experiments, are found and the relative displacement of contacting particles does not depend only on the distance of the centres. For a transversely isotropic aggregate the fluctuation model predicts six independent moduli [20], with no major symmetry of the elastic tensor, *A*_{ijkl} ≠ *A*_{klij}.

## 3. Fluctuation theory

At the end of the section ‘Lattice of multiple point-elements’ [1], Love, in his note B, underlines the need for a theory in which each particle must be in equilibrium under forces exerted upon it by other particles to improve upon the simplest approach. This means that, at the local level, a pair of particles and its neighbours, the configuration is not symmetric in terms of loading and geometry, so equilibrium is needed. This asymmetry can be captured by means of the fluctuations. We review the essential points of the theory proposed in references [9,12] where fluctuations are introduced, and we extend the model to the case of a transversely isotropic aggregate. The kinematics of a pair *A*−*B* is described by
**d**^{(BA)} is the contact vector from the centre of particle *A* to the centre of particle *B*. The increment in contact force is
*A*−*B* have sufficient translational freedom to satisfy force equilibrium, while we assume that the other particles in contact with it translate with the average deformation. Clearly, a more faithful model should also include fluctuations to describe the interaction between the typical pair and its neighbour [13,14]. The result, however, does not lead to a qualitative change in the solution [12] and, therefore, we prefer to work in a simpler context.

For a particle *n* not equal to particle *B*, we assume
*A*. We are aware that this leads to a crude approximation of the equilibrium equations and that including fluctuations in equation (3.3) would lead to a proper formulation for the equilibrium of a typical pair.

The equation of force equilibrium for particle *A* is, then,
*N*^{(A)} is the number of particles in contact with *A*. We determine the solution of the equilibrium equation assuming that *K*_{N} is independent of the contact orientation:
**J** with its average, *Ω*^{(BA)} is the increment of solid angle centred at *M* is the number of pairs in the increment of solid angle. The existence of the contact between particle *A* and *B* provides a symmetry about the plane perpendicular to *A*−*B*. With the conditional averages introduced, equilibrium for particle *B* is also satisfied. A more elaborate solution of the equilibrium of the typical pair *A*−*B* is treated in reference [12] where also the neighbouring particles fluctuate. When employed in equations (3.1) and (3.2), equation (3.9) gives the final expression of the average contact force that particle *B* exerted on particle *A*. The solution for *K*_{N} is constant in the equilibrium equation. However, anisotropy does enter in the problem through equations (3.2) and (2.12). Because *χ*_{1} and *χ*_{2} are given in appendix B. Henceforth, we remove the superscript (*BA*) on the vector **d** for brevity in the exposition. The conditional average of the relative incremental displacement of the centres of particles *A* and *B* can then be derived

Equation (3.12) is the force that particle *B* exerts on particle *A* when both particles, on average, are in equilibrium under the kinematics hypothesis given by equations (3.1) and (3.3). The result of the integration is detailed in appendix C; here, we report the final expression of the elastic tensor
*K*_{N}, equation (2.13), whereas fabric, *β* would differ. The key point in equation (3.14) is that the elastic tensor does not have the major symmetry, *A*_{qlij}≠*A*_{ijql} (or in more detail, *A*_{1133}=*A*_{2233}≠*A*_{3311}=*A*_{3322};*A*_{1122}=*A*_{2211}). Therefore, it is not possible to define a potential for the anisotropic aggregate. This result is achieved by letting particles deviate from the average strain and including anisotropy. If we had *χ*_{2}=0, the contact displacement would have been only a function of the strain applied at a given *β*_{4}=*β*_{5}. Also note that *γ* is a measure of the anisotropy, so *γ*=0 means isotropy.

Another interesting point that arises with the fluctuation model is that the comparison between equation (2.15) and equation (3.14) shows an increment of the number of constants which is consistent with the idea that the simple Cauchy molecular approach leads to a lower number of independent coefficients of the elastic tensor. This is straightforward if we consider the isotropic case, which can be obtained from equation (3.14) by taking *γ*=0, for which only *β*_{2} and *β*_{3} survive [9] (see appendix C). In equation (2.8), there are then two independent coefficients which reconcile the micromechanical model with fluctuations and the continuum approach. Most importantly, the fact that two independent moduli are predicted, instead of the single one of the average strain model, is in agreement with the number of independent moduli of the numerical simulations [9].

As the final step of our activity, we employ numerical simulation to test our finding, that is the lack of major symmetry in the stiffness tensor. This is the only goal in the numerical simulation; we do not make any quantitative comparison with the theoretical model because it has been developed in the simplest form. Yet, despite the simplicity of the model, we anticipate that numerical simulations do indicate the lack of symmetry of the stiffness tensor.

## 4. Numerical simulation

Our goal is to test the lack of symmetry of the stiffness tensor predicted by the analytical model, through numerical simulations. We carry out a numerical simulation [21] to generate a transversely isotropic aggregate made of 10 000 frictionless, elastic, spheres with radius *R*=0.1×10^{−3} m. The shear modulus of the material of the spheres is *μ*=2.9×10^{10} Pa, and the Poisson ratio is *ν*=0.2. Particles interact via normal forces that follow the Hertz law. The initial state is obtained in the manner described in reference [22]. An initial random aggregate of frictionless spheres without contacts is homogeneously and isotropically contracted, bringing the spheres into contact, until a pressure *p* of 976×10^{3} Pa is reached. In this state, the coordination number *k*=6.33 and the solid volume fraction *v*=0.6423. Then, the material is strained along one direction, say *y*_{3}≡**h,** keeping the pressure constant [23]. A triaxial compression in a pure frictionless system is quite complicated, because the system tends to be unstable; so we introduce a friction coefficient equal to 1×10^{−3} to reach equilibriated states as described in appendix D. A careful process of straining is applied: the axial, incremental strain is of the order of 10^{−7} and after each increment, the system is relaxed towards the new equilibrium state. During the straining, we apply *E*_{33}<0, with resulting positive strains *E*_{11}=*E*_{22}, as the pressure is kept constant. We focus our attention on regime of deformations where the shear strain *γ* is small compared with the volume, *ε*_{0}, associated with the pressure. In the initial, isotropic state, we measure
*A*_{1212}=34 MPa, *A*_{1313}=27 MPa, *A*_{2323}=26 MPa.

Some considerations follow: the initial state is almost isotropic, but the response is not exactly the same in all three directions (the *y*_{3}-direction seems little bit stiffer); the moduli are evaluated within a range of deformation 1×10^{−6} and 9×10^{−6}; the presence of a small friction coefficient, 1×10^{−3}, does not affect the solution with respect to the pure frictionless condition.

Computer simulations [9,24] suggest that the average strain theory predicts reasonably well the bulk modulus

In figure 1, we plot the deviatoric strain *γ* normalized by *ε*_{0} versus the deviatoric stress *p*. In the initial state, the dimensions of the box *L*_{1}, *L*_{2}, *L*_{3} are identical, whereas *σ*_{3} is slightly bigger than *σ*_{1} and *σ*_{2}. Therefore, at zero deviatoric strain, we have a non-zero deviatoric stress. During the axial compression, there is a small, negligible, variation of both the coordination number and the fabric. The anisotropy is mainly owing to the elastic deformation of the grains and it ensures an induced anisotropy in the aggregate. In figure 1, we also show the stressed states where probes are applied, apart from the initial state. The incremental response is achieved by applying an incremental unloading with respect to the previous forward loading along the axial strain. Details of the simulations are given in appendix D. So, for example, as we compress along *y*_{3}≡**h** with *E*_{33}<0 during the axial strain, in order to evaluate *E*_{33}<0, the pressure is kept constant and we have *E*_{11}=*E*_{22}>0. Therefore, the incremental response

We report the data associated with probe n.1 (figure 1):

There is a clear evidence of the difference between *A*_{1133} and *A*_{3311}, *A*_{2233} and *A*_{3322}.

The theoretical prediction of the incremental response can be derived from equation (3.14). In particular, in the initial state, where *k*=6.32, *v*=0.6423, *p*=976 kPa, we have
*γ*/*ϵ*_{0}=0.012, *p*=976 kPa, *k*=6.32)
*γ*/*ϵ*_{0}=0.027, *p*=976 kPa, *k*=6.33)

The theory predicts a rather small difference between the moduli that we identify in two main reasons. The first is that in the analytical model the difference (*ρkχ*_{2}*ϵ*^{1/2}/5)×*γ*/*ϵ* is proportional to *γ*/*ε* which is small in the regime of deformation of our interest. A frictionless aggregate cannot be strained enough to induce a strong anisotropy before failure occurs (a stronger difference is possible in the frictional case [25]). The second reason is due to the simple model adopted. The difference in the moduli is related to the fluctuations. Several works [9,13,14] have pointed out the need to include more degrees of freedom in the kinematics of contacting particles in order to capture a more faithful response of the aggregate, with a correlation length that increases as we approach the isostatic condition of the aggregate, *k*=6 for a frictionless case [26]. So, because we adopt a simple pair fluctuation model and our aggregate is close to the isostatic condition, we expect to capture a qualitative feature rather than a quantitative prediction. On the other hand, if we employed a model with more degrees of freedom, we would have dealt with a far more complicated theory at the expense of the simplicity here developed. We also note that while the theory predicts *A*_{3311} greater than *A*_{1133}, numerical simulations show the opposite.

An energetic issue arises with the lack of major symmetry of the stiffness tensor. While we are in the process of a complete analysis about this point, here we suggest that because of the working of the fluctuations there is no potential energy for the average strains.

Finally, we have also carried out numerical simulations using a linear contact model. There, anisotropy is given by the fabric [17]. We recover similar results to the Hertizan contact model, and confirm the lack of major symmetry of the stiffness tensor. We underline that it is crucial to restrict the attention in the pre-peak region where homogeneous deformations are likely to occur and where the present theory holds.

## 5. Conclusion

A micromechanical analysis for a random, aggregate of frictionless, elastic particles has been employed. The model differs from that based upon an affine deformation because of the presence of a deviation that ensures pair particles equilibrium. In the transversely isotropic case, we derive an elastic tensor with six constants, instead of the expected five, and without major symmetry. We achieve this result because of the simultaneous presence of anisotropy and fluctuations in particles kinematics. Numerical simulations confirm our finding.

## Data accessibility

This work does not have any experimental data. We generate numerical data by means of an update version of Trubal code, distinct element method.

## Author' contributions

L.L.R., G.R. and A.S. employed the theory. L.L.R. and L.O. carried out the numerical simulations. All authors approved the final version of the manuscript.

## Competing interests

We have no competing interests.

## Funding

L.L.R. is grateful to GNFM ‘Gruppo Nazionale della Fisica Matematica’ (Italy) for the support. He is also grateful to ONR Global through grant no. N62909 13 1 V152. L.O. and L.L.R thank National Science Foundation grant no. PHY11-25915 for support of their visit at the Kavli Institute of Theoretical Physics in Santa Barbara, CA, USA.

## Acknowledgements

We grateful to the referee for his suggestions and corrections on numerical simulations.

## Appendix A

When we apply the divergence theorem to tensor products on the unit sphere, we obtain the following identities:

The elastic tensor is

## Appendix B

By definition, the continuous version of the average of tensor **J** is
*B*. We follow [9] and we have

## Appendix C

The incremental stress is

So, we have
**X**,**Y** and **Z** defined in appendix A. We obtain
*γ*=0,

## Appendix D

We refer to the anisotropic state, probe n.2, and we report numerical data to show in details our results. Before we apply any increments to determine the moduli, we stabilize the anisotropic state. We introduce a quantity proportional to the average kinetic energy
**v** is the particle velocity. Then, we consider a dimensionless quantity
*t* is the numerical time step and the average displacement in the initial state is
*R* is stable over a long numerical time in which no strain is applied to the box.

In figures 3 and 4, we show the evolution of the pressure and the coordination number in the same anisotropic state over a long time of relaxation, before any increments is applied. There is again evidence of a clear stable numerical state.

The incremental response is determined by applying an affine motion followed by a long relaxation time. In figures 5 and 6, we refer, respectively, to the evolution of the coordination number and the pressure when a deformation 6.8×10^{−6} is applied.

In figure 7, we plot the evolution of the effective moduli *A*_{1122} and *A*_{2211} under different relaxation time for different applied strain (figure 8).

In figure 9, we plot the evolution of the effective moduli *A*_{1133} and *A*_{3311} under different relaxation time for different applied strain (figure 10).

In figure 11, we show, as example, the evolution of *A*_{1133} with the applied strain, from 10^{−7} to about 10^{−5}.

- Received January 8, 2015.
- Accepted October 1, 2015.

- © 2015 The Author(s)