Internal states of model isotropic granular packings.II. Compression and pressure cycles.

Internal states of model isotropic granular packings. II. Compression and pressure cycles.


This is the second paper of a series of three investigating, by numerical means, the geometric and mechanical properties of spherical bead packings under isotropic stresses. We study the effects of varying the applied pressure (from 1 or 10 kPa up to 100 MPa in the case of glass beads) on several types of configurations assembled by different procedures, as reported in the preceding paper Agnolin and Roux (2007a). As functions of , we monitor changes in solid fraction , coordination number , proportion of rattlers (grains carrying no force) , the distribution of normal forces, the level of friction mobilization, and the distribution of near neighbor distances. Assuming that the contact law does not involve material plasticity or damage, is found to vary very nearly reversibly with in an isotropic compression cycle, but all other quantities, due to the frictional hysteresis of contact forces, change irreversibly. In particular, initial low states with high coordination numbers lose many contacts in a compression cycle, and end up with values of and close to those of the most poorly coordinated initial configurations. Proportional load variations which do not entail notable configuration changes can therefore nevertheless significantly affect contact networks of granular packings in quasistatic conditions.

45.70.-n, 83.80.Fg, 46.65.+g, 62.20.Fe

I Introduction

The mechanical properties of solidlike granular packings are traditionally studied, at the macroscopic level, in engineering fields such as soil mechanics Wood (1990); Biarez and Hicher (1993); Mitchell (1993); Herrmann et al. (1998), and are currently being investigated, with some attention to the grain scale and micromechanical origins of macroscopic behaviors, in condensed matter physics and material science communities Herrmann et al. (1998); Hinrichsen and Wolf (2004); García Rojo et al. (2005).

The present paper, the second of a series of three, investigates, by numerical simulations, the mechanical and microstructural response of a model material, the packing of identical spherical beads, to pressure intensity variations. It refers a lot to the results of the previous, companion paper Agnolin and Roux (2007a), but may be read independently.

Although molecular dynamics (or “discrete element”) approaches have repeatedly been applied to sphere packings Makse et al. (1999); Thornton (2000); Makse et al. (2000); Silbert et al. (2002a, b); Suiker and Fleck (2004), many important questions related to the microscopic origins of their macroscopic mechanical behavior in the quasistatic regime have not been fully explored yet. One such issue is the influence of the initial state, which is determined by the assembling process. In the first paper of the present series Agnolin and Roux (2007a) (hereafter referred to as paper I), the results of several packing preparation methods, all producing ideally isotropic states, are compared. Direct compressions of granular gases produce states that do not depend on dynamical parameters if the compression is slow enough. Their solid fraction and coordination number (evaluated on excluding the rattlers, a proportion of grains which do not carry any force) are decreasing functions of the friction coefficient , from and for , in which case the random close packing state (RCP) is obtained, down to and for . In paper I Agnolin and Roux (2007a) we accurately checked the uniqueness of the RCP, on confronting our own numerical results with those of several recent publications, in which different numerical procedures were implemented Donev et al. (2005); O’Hern et al. (2003). In the presence of intergranular friction, however, quite different packing states might be prepared. First, it is of course possible, in a simulation, to increase the friction coefficient once the packing is equilibrated under some pressure; such a numerical procedure can be regarded as a model for an assembling process in the presence of a lubricant within intergranular gaps in the laboratory. Ideally, whatever the value of the friction coefficient used to model the quasistatic mechanical properties of the material, it is possible to assemble the sample with (thus assuming ideal, perfect lubrication in the fabrication stage) and hence with the RCP density and coordination number. Once the grains are packed and form a solid material, contacts between grains can then be attributed the final, finite friction coefficient used in quasistatic modelling. Experimentally, it is of course well known that given granular materials can be packed with varying densities. A common method to make them denser, other then lubricating the contacts in the assembling stage, is the application of vibrations or “taps”. A numerical idealized vibration procedure, apt to prepare dense samples with little computation time, was defined in paper I. Surprisingly, although it produces isotropic states with densities close to the RCP value, their coordination numbers are as low as in the loosest states assembled by direct compression. The small geometric differences between configurations with the same solid fraction but very different coordination numbers is still not accessible to tomographic observation techniques Agnolin and Roux (2007a). Only mechanical properties can thus be confronted to experimental results, to determine whether or in which conditions the investigated numerical systems are close to experimental reality.

Before studying elastic properties in paper III of the present series Agnolin and Roux (2007b), one should first investigate the effect of an isotropic compression. The application of a large enough confining pressure, usually at least a few tens of kPa (with rare exceptions Tatsuoka et al. (1986); Reydellet and Clément (2001)), is necessary before the macroscopic mechanical behavior of solidlike granular packings is tested Vermeer (1998); Gudehus et al. (1984); Wood (1990), and characteristic quantities such as dilatancy and internal friction angle are measured. Experimental data on elastic moduli Thomann and Hryciw (1990); Shibuya et al. (1992); Jia et al. (1999); Jia and Mills (2001); Kuwano and Jardine (2002); Geoffroy et al. (2003); Sharifipour et al. (2004); Agnolin et al. (2005) are also extremely scarce below that range. Most relevant laboratory sample histories to be understood in order to relate the macroscopic response to internal variables and micromechanics involve an assembling stage, and then a compression stage, which is often isotropic or oedometric. It is therefore necessary to assess the influence of pressure changes on the initial states.

In addition, the material behavior under varying isotropic stress is interesting per se. The behavior of sands is traditionally regarded Wood (1990); Gudehus et al. (1984); Biarez and Hicher (1993) as elastoplastic under isotropic loading, with pressure cycles entailing irreversible density increases. Such effects are nevertheless considerably smaller than in cohesive materials such as clays Wood (1990); Biarez and Hicher (1993); Mitchell (1993), or powders Castellanos (2005). It is worth investigating such behavior in model sphere packings by numerical means.

Ii Model material, micromechanical parameters

ii.1 Contact model

We briefly recall here the model material and the contact laws, which are described in paper I with more details. Equal-sized spherical beads of diameter (whose value, as we ignore gravity, will prove irrelevant), interact in their contacts by point forces of elastic, frictional and viscous origins. The Hertz law relates the normal elastic force to the normal deflection (approach of sphere centers closer than ) as :


with the notation , being the Young modulus of the beads, and the Poisson ratio. The Hertz law introduces a normal stiffness that depends on or on .

Tangential elasticity and friction are described with a simplified form of the Cattaneo-Mindlin-Deresiewicz results Johnson (1985), in which the tangential stiffness , relating the tangential elastic force increment to the relative tangential elastic displacement in the contact, is proportional to :


The Coulomb condition with friction coefficient requires to be projected back onto the circle of radius in the tangential plane whenever the increment given by Eqn. (2) would cause its magnitude to exceed this limit. In order to avoid unphysical increases of elastic energy, is scaled down in proportion with when the elastic normal force decreases, as indicated in paper I and advocated in Elata and Berryman (1996). Tangential contact forces also move with the particles in contact, so that the condition of objectivity is satisfied (see paper I and ref. Kuhn and Chang (2006)).

A viscous term opposing normal relative displacements reads (positive normal forces are conventionally repulsive):


with a damping coefficient depending on elastic normal deflection (or on elastic repulsive force ), such that its value is a fixed fraction of the critical damping coefficient of the normal (linear) spring of stiffness joining two beads of mass :


We do not introduce any tangential viscous force, and impose the Coulomb inequality to elastic force components only. The main justification of such a term is computational convenience (to accelerate the approach of equilibrium states), and we could check that its value did not affect the statistical results on the configurations of the packings.

The present numerical study was carried out with the elastic parameters and that are suitable for glass beads, and the friction coefficient is set to .

ii.2 Stress control

The numerical results presented below were obtained on samples of beads, enclosed in a cubic or parallelipipedic cell with periodic boundary conditions. The sizes of the cell are denoted as , parallel to coordinate axes (). ’s vary simultaneously with the grain positions and orientations until mechanical equilibrium of all particles with the prescribed values of all three diagonal components of the Cauchy stress tensor, , is obtained. One then has :


Here is the sample volume, ’s are the coordinates of vector joining the center of bead to the one of its contacting neighbor (with the nearest image convention of periodic cells) and ’s are those of the corresponding contact force. This force is actually exerted by onto , so that the convention used is that tensile stresses are negative. Velocities of grain centers comprise, in addition to a periodic field, an affine term corresponding to the global strain rate. Equations of motion for dimensions are written in addition to the ordinary equations for the dynamics of a collection of solid objects, and they drive the system towards an equilibrium state in which condition (5) is obeyed.

In the present study we always impose isotropic stresses, i.e. hydrostatic pressures : for , 2, 3.

ii.3 Dimensionless parameters

In addition to include friction coefficient and viscous dissipation parameter , the important dimensionless control parameters for sphere packings under given pressure are the reduced stiffness and the inertia parameter . is chosen such that the typical contact deflection is proportional to ,


a correspondance which can be made accurate thanks to the relation


between pressure and the average normal force in the contacts. (7) is exact provided in all contacts and intercenter distances are taken equal to the diameter . Here denotes the cordination number, equal to , with the total number of force-carrying contacts in the packing. Rattlers, in proportion , have no such contact. We refer to te force-carrying network - the packing devoid of its rattlers – as the backbone, and to , which simply relates to as r, as the backbone coordination number. Brackets denoting averages over all force-carrying contacts, one has

The limit of rigid grains is approached as .

can be used to determine whether the material within the grains is likely to be imposed stresses beyond its elastic limit. The maximum pressure, at the center of a Hertzian contact between spheres of diameter , carrying a normal force , is Johnson (1985)

Under pressure , corresponding to by (6), when the average normal force in contacts is , one can deduce from (7)


Likewise, the maximum shear stress , which is reached inside the grains near the contact region will be Johnson (1985) (for )


Eqns. 8 and 9 show that very high stress levels, up to a non-negligible fraction of elastic modulus are reached if is not large enough. With our choice of material parameters for glass beads, we get for P=10 MPa and for P=100 MPa, while the numerical prefactor is only slightly lower than 1 () if (a typical value) in (8). Such high stresses are very likely to entail particle breakage or plastic strains (according to the materials the grains are made of).

In our simulations we set our lowest pressure level for the simulation of glass beads to 1 kPa or 10 kPa, corresponding to and with the elastic properties of glass. This enables us to explore the entire experimental pressure range, and to approach the large limit too. Up to the maximum pressure value 100 MPa, we assume elastic contact behavior, but one should be careful on comparing the numerical results in the higher pressure states ( MPa) to experimental ones.

Dynamical effects are assessed on comparing the strain rate to intrinsic inertial times, such as the time needed for a particle of mass , initially at rest, accelerated by a typical force , to move on a distance . This leads to the definition of a dimensionless inertia parameter :


The quasistatic limit can be defined as . is a convenient parameter to describe internal states and write down constitutive laws for granular materials in dense shear flow MIDI (2004); da Cruz et al. (2005); Pouliquen et al. (2005); Jop et al. (2005, 2006).

ii.4 Initial states

Procedure (%) Z(2)
A 0 0
B ()
C (vibration)
Table 1: Isotropic states ( for A and C, for B and D) for different assembling procedures.

The present paper is devoted to the study of the influence of quasistatic pressure changes to granular packings assembled by different means, as described in paper I Agnolin and Roux (2007a). Four different states were prepared under low pressure, and some of their basic characteristics are recalled in table 1. Such state variables are monitored in the following as a function of pressure in isotropic compression or pressure cycles. In addition to solid fraction , proportion of rattlers , backbone (or force-carying structure) coordination number , Table 1 provides some global information on force distributions. is characteristic of the width of the distribution of normal forces:


and are the average levels of friction mobilization (i.e., ) for contacts carrying normal forces, respectively, larger and smaller than the average .

In paper I we also recorded other geometric data, in particular pair correlation functions and distributions of near neighbor gaps . The latter can be expressed as gap-dependent coordination numbers, defining as the average number of neighboring beads around a central one, separated by an interstice smaller than . thus coincides with the contact coordination number. Due to the rattlers, the proportion of which –see table 1– can exceed 10% of the total number of grains, such geometric data are however somewhat ambiguously defined: the positions of the rattlers are not fixed by the rigid backbone. Thus one may define , on using the arbitrary positions obtained at the end of the simulation, when the packing first equilibrates within the prescribed numerical tolerance. One then has (recall counts only force-carrying contacts) if the equilibrium state is accurately computed, because there are very few contacts bearing a normal force below tolerance. In an attempt to define more intrinsic geometric data, we defined in paper I Agnolin and Roux (2007a) as the gap-dependent coordination number in the configuration obtained once all rattlers are pushed against the backbone, in random directions. In their new position, the rattlers now have three contacts with the backbone (except in the rare case when inter-rattler contacts are obtained). It was argued in paper I that the resulting structure was likely to resemble, to some extent, granular assemblies under gravity, when the weight of the grains is very small in comparison to the local stress. can be regarded as a geometric definition of a contact coordination number (it is, in general, slightly larger than ).

Iii Numerical results

We first specify the numerical compression procedure in paragraph III.1, then describe the effects of an isotropic compression and a pressure cycle in terms of global variables (Section III.2) as well as local geometry (Section III.3). We then test the simplest prediction scheme for the evolution of coordination number, that of homogeneous strain at the microscopic level, in Section III.4.

iii.1 Numerical procedure

The results presented below pertain to equilibrium configurations at variable isotropic pressure , obtained by a stepwise compression (respectively: decompression) process in which , within the controlled stress scheme described in Section II.2, is increased (respectively: decreased) by a factor . In each pressure step a condition of slow enough strain rate was enforced, so that the inertia parameter, as defined by (10) with the currently imposed pressure level, was kept below a maximum value: for compression, for decompression. Such values were chosen to ensure independence of the results on dynamical parameters and . It was observed that a decompression process requested more care, due to its greater instability. Whereas a compression of the sample beyond its equilibrium density will be strongly opposed, at growing , by elastic forces in the network, too large an expansion, as decreases, might cause the contact network to break apart, resulting in a dynamical process similar to assembling a granular gas, when the externally applied pressure finally drives the system back to a denser equilibrium configuration. Such events might entail a significant remoulding of the contact network and large departures from equilibrium conditions. This should of course be avoided in a procedure designed to model a quasistatic evolution, as close as possible to the limit of small strain rates.

Configurations are deemed equilibrated when, defining as a small tolerance on forces and as a small tolerance on energies, the four following conditions are simultaneously satisfied :

  • each coordinate of the total force on each grain is smaller than ;

  • each coordinate of the total moment on each grain is smaller than ;

  • all stresses have their prescribed values with a relative error smaller than :

  • the kinetic energy per grain is smaller than .

To distinguish between the backbone and the rattlers, the same method is applied as presented in paper I Agnolin and Roux (2007a).

Such procedures were applied to samples A to D below, with ranging from its smallest value 1 kPa (for B and D, corresponding too ), or 10 kPa (for A and C, corresponding to ), up to 100 MPa (), and then back to its initial low value. Letters A, B, C, D will hereafter denote pressure-dependent configuration series. Although initial states A and B were assembled with coefficients of friction lower than the chosen value , we study quasistatic compressions with for all sample series. We regard the smaller friction levels applied to configurations A and B in the assembling stage as models for lubricated grains, and assume that the lubricant ceases to operate once solid particles finally touch one another, as in equilibrated packings and during quasistatic compression tests. As a reference for comparisons with other states, and because it was studied in the literature Makse et al. (2000, 2004), we also prepared another configuration series we denote as A0, obtained from the initial A state on compressing a frictionless system (thus series A and A0 share the same initial low-pressure state, but differ as soon as is altered).

All results are averaged over 5 samples of n=4000 beads, and error bars correspond to one standard deviation.

iii.2 Evolution of global state variables

Figs. 1 and b display the evolution of solid fraction , backbone coordination number , and rattler fraction in sample series A, B C and D in the pressure cycle.

Figure 1: (Color online) Evolution of packing fraction as a function of pressure in glass bead packings (bottom axis), or dimensionless stiffness parameter (top axis), in (from top to bottom) states A (red crosses, continuous line), C (black square dots, continuous line) B (blue asterisks, dotted line) and D (green open squares, dotted line).

Fig. 1 shows that the solid fraction change with pressure is almost perfectly reversible: the data points corresponding to the compression and decompression parts of the pressure cycle are almost indistinguishable. More precisely, once the pressure had returned to its lowest value in samples A to C, the packing fraction was observed to have changed by very small amounts, below . The loosest state, D, undergoes a slight compaction. Yet, this effect apparently decreased as the maximum prescribed value for parameter was changed from to upon unloading (the reported results corresponding to this latter value). Our model material thus differs from sands, which are reported to respond to such cycles with notable irreversible density increases Biarez and Hicher (1993); Wood (1990). It should be noted, though, that we are using a contact model without plasticity or particle damage, which, as argued on evaluating, in Sec. II.3, the maximum pressure and shear stress in the grains near contact points with Eqns. (8) and (9), is quite unrealistic for the highest pressure levels simulated. Stress concentrations in contacts between angular particles like sand grains, with corners or asperities Johnson (1985); Goddard (1990), are more severe than between smooth objects and should enhance the effects of anelastic material behavior within the grains. The smallness of irreversible compaction in our simulations suggests that such macroscopic behavior, in sands, originates in contact mechanics rather than in collective effects.

The reversibility of the response to the pressure cycle is however only apparent, as the coordination number does not return to its initial value.

As expected, increases under a growing confining pressure (Fig. a): as the particles pack more closely in a smaller volume, near neighbors come into contact. reaches about 7.3 at the highest pressure in the densest samples, A and C. Correlatively, an increasing number of rattlers get trapped as their free volume shrinks, and are recruited by the force-carrying network. The initially large fraction of rattlers in states C and D () steadily decreases as grows( Fig. b) and has virtually disappeared at  MPa.

(a) versus or
(b) versus or
Figure 2: (Color online) Backbone coordination number (a) and proportion of rattlers (b) as functions of or , same symbols as on Fig. 1.

The evolution of coordination numbers on unloading is more surprising. While low coordination states C and D exhibit a very limited hysteresis effect and eventually retrieve their initial, low values (about ), with a slightly lower rattler fraction, samples of types A and B, in which was initially high, lose contacts as a result of the pressure cycle and end up with values below 5 (about 4.8 for A, and 4.5 for B), closer to C and D ones than to where they started, with a substantial rise in the population of rattlers. (Let us recall that samples A and B are regarded in the study of quasistatic compression as made of frictional beads with , like the others). The behavior of (frictionless) samples A0 is of course different, for they cannot be stable at low pressure below  Roux (2000). Fig. 3 compares the evolutions of in A and A0 series, and shows that is very nearly reversible in the A0 series. The unloading curves in A states starting at lower pressures, 3.16 Mpa and 1 MPa instead of 100 MPa, also shown on Fig. 3, witness a lower, but significant decrease of from its initial value at the end of the cycle.

Figure 3: versus or in pressure cycle in series A (crosses) and A0 (dots), showing reversibility for A0. Shorter cycles (up to 0.316 MPa and 1 MPa) than the one of Fig. b are also shown for A.

The shape of the force distribution and the mobilization of friction also change with , as shown by the evolution of parameters , , on Fig. 4.

Figure 4: (Color online) From bottom to top: , and versus or in compression cycle. Symbols as on Fig. 1 for states A, C, and D. Series A0 represented with (red) dots joined by dotted line for . Hysteresis loops for first decrease, then increase back on unloading and go through a maximum (except for A0, in which cas it is nearly constant). and behave in a similar way, with the special circumstance that their initial values are equal to zero in A states (assembled without friction).

As a general rule, the width of the force distribution correlates with the level of force indeterminacy, relatively to the number of degrees of freedom. Contact elasticity tends to share forces rather evenly, because contact force values should minimize the intergranular elastic energy, subject to the constraint that they balance the applied pressure (this elastic energy as a function of forces is written further below in connection with a discussion of irreversibility in pressure cycles, and the minimization property is exploited in paper III Agnolin and Roux (2007b) to estimate bulk moduli). More precisely, the increments of forces due to pressure increases will tend to reduce the width of the distribution, the faster the less constrained the minimization, i.e. the larger the degree of force indeterminacy. Thus in configurations A, the large coordination number enables a quick narrowing of the distribution under growing pressure. In states C, the same tendency is present, but the evolution is much slower, as there are less possibilities to distribute forces in a more tenuous network while maintaining equilibrium. However, C samples gain contacts faster than D ones (Fig. a), for which the narrowing effect is even slower. Finally, the extreme case is the situation of isostaticity, as in the A0 series, in which the distribution of forces is geometrically determined in the rigid limit of . As, furthermore, the increase of with is not very fast in that case, since is already large from the beginning, the shape of the distribution remains nearly constant. A few normal force probability distribution functions at different pressure levels are shown on Fig. 5.

Figure 5: From bottom to top, evolution of normalized force distributions , with , with growing pressure in samples A, C, D, A0. value in kPa are (except for D: ), , , and . All four distributions tend to narrow as grows, but at very different rates.

The evolution of force values and friction mobilization on unloading is more complicated: all three parameters shown on Fig. 4 first increase, then go through a maximum and end up, at the initial pressure value, with a value comparable with the initial one (except for friction mobilization parameters and in A systems, because they started at zero). In a granular sample controlled in displacements or strains, rather than stresses, large self balanced forces can in some situations remain when the external load that created them is removed, the simplest example being that of one particle wedged in a corner McNamara et al. (2005); Halsey and Ertaş (1999). Our observations indicate that such a phenomenon does not take place in a situation of controlled stress state: all forces are of the order of the average force, which is related to the current pressure by (7), even though contacts have carried forces that were larger by orders of magnitude in the past. This suggests that the set of admissible contact forces, restricted to the intersection of an -dimensional affine space (due to equilibrium relations) with a cone (due to Coulomb inequalities) is bounded. Yet during unloading many more sliding contacts are observed than at growing pressure, due to the effects of decreasing normal force components, and the level of friction mobilization is higher (Fig. 4). Meanwhile, the distribution of normal forces gets wider. The global influence of the past loading, with contacts previously carrying larger forces, enhances force heterogeneities. A related quantity is the elastic energy stored in the contacts. The total elastic energy per grain reads (from Eqns. (1) and (2))

Once adimensionalized by , we denote it as . On exploiting Eqn. (7) it is conveniently expressed as:


In (12), , related to force moments, is close to , which can be defined on replacing exponent 2 by in Eqn. (11), with the following slight modification. With defined in (2) as the constant ratio of tangential to normal stiffnesses, and with the notation for the ratio in a contact, let us define


thus depends on the force distribution and also on friction mobilization, although for its relative difference with is small (of the order of , with as plotted on Fig. 4). The energy per particle, , scales as , which is expected since this is proportional to for the typical normal contact deflection. is larger for low coordination numbers (weaker networks), and larger force disorder (higher ). (It should be recalled that we use pressure, rather than strain, as the control parameter, hence a larger elastic energy for softer materials). Thus in A configurations, is larger, for given , on decompressing, another manifestation of the irreversibility of the cycle. If we assume that the curve is quasistatically followed up to the maximum pressure, and then exactly retraced back on decompressing, this leads to a paradox, as some elastic energy appears to be gained at no expense. Thus one has to account for very small irreversible density changes, for energetic consistency. Such changes in , between the growing and decreasing pressure parts of the cycle are shown on Fig. 6.

Figure 6: (Color online) Increment of packing fraction gained between the two states of equal pressure, reached at growing and at decreasing , in states A (red, crosses) and C (black, square dots). Note the scale of density changes ( of order ).

In the case of A configurations, one even observes a slight decompaction on decreasing back to its lowest, initial value. Although surprising, this phenomenon should be expected in the rigid limit or , because as explained in paper I, the initial A configuration, which was assembled without friction, is a local maximum of subject to impenetrability constraints. Another conclusion of paper I Agnolin and Roux (2007a) is that the only way to increase density in such a sample is to produce, by enduring agitation or repeated shakes, notable traces of crystalline order. This should not happen in a slow, quasistatic compression experiment with only one pressure cycle. To check for energetic consistency, one may note on Fig. 6 however that the change of is positive at high pressure. The total energy fed into the system in the cycle is


the integral running over the whole pressure interval of the compression cycle. Consequently (see Fig. 6) the contribution of the irreversible increase of is largely dominant, because it is integrated over a much wider pressure interval. The small changes in density between the compression and the decompression curves at the same pressure values are large enough to explain the change in elastic energy, and that of potential energy as well when the cycle ends up decreasing the density (which happens for A samples).

iii.3 Pair correlations and near neighbor distances

The smallness or absence of irreversible compaction in the pressure cycle implies that the samples do not avoid contact deflections by finding denser packing arrangements. Thus interparticle correlation patterns should witness favored near neighbor distances which typically scale like .

(a) versus in C configurations at different .
(b) versus in C configurations at different .
Figure 7: Pair correlation functions at P=10, , , ,  kPa in configurations C at growing pressure, without (top), and with (bottom) rescaling distance as .

This is shown for C configurations on Fig. b: on rescaling the distance axis, using coordinate with the initial low pressure solid fraction, the different curves are superimposed. In agreement with the observations made in paper I Agnolin and Roux (2007a), where the relationships between pair correlation functions and contact networks were discussed, a closer look on such correlations will reveal differences in the details of the peaks associated with changes in the coordination number with . Figs. 8 and 9 respectively show functions and at growing values, using the corresponding change of scale for interstice , .

Figure 8: (Color online) Gap-dependent neigbor coordination number versus rescaled interstice at different (same as on Fig. b) in states A (red), C (black) and D ( green).
Figure 9: (Color online) Same as Fig. 8 for definition of the gap-dependent neigbor coordination number.

Those data suggest that the homogeneous shrinking of distances implied by the rescaling of abscissae on the graphs of Figs. b, 8 and 9 is an approximation with some discrepancies at small intergranular distances. Curves corresponding to pressures other than the lowest one on Figs. 8 and 9 start at distance and the corresponding values of on the curve for the lowest pressure value are the predictions for the coordination number on assuming homogeneous shrinking strains. Differences therefore show that such predictions, albeit reasonable, are not exact. In particular, the gradual capture of rattlers by the force-carrying network as grows (see Fig. b) cannot be adequately described by the homogeneous shrinking assumption: the rattlers will not start carrying forces when one interstice with a backbone grain is closed. The use of definition should in principle improve this kind of prediction: once positioned against the backbone (with 3 contacts), the rattlers are much more likely to create new contacts bearing nonzero forces when they touch new neighbors. Yet, the improvement of curve superpositions on Fig. 9 compared to Fig. 8 is marginal. This suggests that the inaccuracy of the prediction of coordination numbers is not only due to the capture of rattlers by the growing backbone, but also stems from the failure of the assumption of homogeneous shrinking.

iii.4 Can one predict the changes in coordination number ?

The results of the prediction of the coordination number, assuming all distances uniformy shrink, are shown on Fig. 10 for systems A and C under growing pressure. The agreement is very good in state A (except at high pressure, where is slightly underestimated), and fair in state C. For C configurations, the prediction was done separately for both and , showing a somewhat better accuracy at low pressure in the second case. Unfortunately, the mechanically important coordination number is rather than .

Figure 10: (Color online) Predictions for in samples A and C, and for in samples C, based on the homogeneous shrinking assumption.

To evaluate as a function of , one needs to account for two phenomena: the increase of the elastic normal deflection in the contacts that already existed at the lowest pressure, and the creation of new contacts due to the closing of open interstices. Both effects are evaluated with the assumption of homogeneous rescaling of all distances according to the density change, respectively exploiting the previous measurements of the distribution of sphere overlaps (related to that of normal forces), and of the function (with no significant difference in accuracy on using or ). The predicted values of , although not very accurate for small changes of at low pressures, globally capture the marked growing trend above 1 MPa. The predictions of density increases are compared with the simulation results on Fig. 11, showing good agreement (with a slight underestimation at high pressure). The prediction of is understandably more accurate than that of the coordination number, because it is not very sensitive, at first, to errors in the estimation of the density of newly created contacts, which initially carry very small forces.

Figure 11: (Color online) versus or in samples A (red) and C (black). Dots: measurements. Dotted lines: predictions, based on the homogeneous shrinking assumption from the initial state of lowest pressure.

One may also attempt to predict the decrease of coordination number in the decompression part of the pressure cycle. Such a prediction is based on the distribution of particle overlaps (or contact deflections), rather than near neighbor distances. The relevant information is therefore the normal force histogram for the highest pressure level, as shown, e.g. on Fig. 5. However, this is a rather crude approximation, which leads to large errors for the coordination number variation with density, as shown on Fig. 12, and very poor predictions indeed for the coordination number relationship to the decreasing pressure, as apparent on Fig. 13.

Figure 12: (Color online) Coordination number versus at decreasing in samples A (red) and C (black). Dots: measurements. Dotted lines: predictions, based on the homogeneous expansion assumption from the initial state of highest pressure.

Such an assumption of homogeneous expansion proves in particular unable to provide a correct estimate of the properties at low density or pressure, as it ignores the requirement of mechanical rigidity.

Figure 13: (Color online) Coordination number versus decreasing (or ) in samples A (red) and C (black). Dots: measurements. Dotted lines: predictions, based on the homogeneous expansion assumption from the initial state of highest pressure.

We are not aware of a simple prediction scheme that would be able to provide a reasonably accurate description the reduction of coordination number in the A state on reducing the confining pressure.

Iv Discussion

The effect of a compression on the four series of isotropic packings we have been studying can be broadly summarized as the closing of additional contacts and the gradual reduction of the characteristic disorder of granular systems, as witnessed by the narrowing of the force distribution (Figs. 4 and 5). Geometric changes conform to the homogeneous shrinking assumption on large scale, and the resulting predictions for the near-neighbor distances and the coordination numbers are reasonable, if not very accurate, approximations (Figs. 10 and 11), even though they cannot correctly account for the recruitment of rattlers (Fig. b) by the growing backbone. It proves difficult to accurately estimate small increases, to which, as will be studied in Agnolin and Roux (2007b) (paper III), shear moduli of poorly coordinated packings are especially sensitive. The changes in the forces and the mobilization of friction are not appropriately described by such a simple model. On assessing the performance of the homogeneous shrinking approximation, one thus retrieves the classification of length scales introduced in paper I (Agnolin and Roux, 2007a, Section IV.E.2). Global changes on scales above about appear to abide by the homogeneous strain assumption, hence the superposition of pair correlation functions on Fig. b. Pair correlations between neighbors at smaller distances (or details of the peaks of ) are only approximately predicted on rescaling all distances by the same factor (as appears on Figs. 8 and 9). And small distances of the order of (contact deflections related to forces) do not abide by this homogeneity of strain. Otherwise, on rescaling coordinates by a factor , where , one would replace any contact deflection by , which for would result in a much stronger narrowing of the force distribution than the one observed. This assumption of homogeneous strain (or affine displacements) will be further tested on dealing with elastic moduli in paper III Agnolin and Roux (2007b).

The effects of a pressure reduction are more surprising. Although the evolution of solid fractions departs very little from reversibility (Figs. 1 and 6), large initial coordination numbers in configurations A and B do not survive a pressure cycle (see Figs. a and 3). Such effects are not predicted by the simple assumption of homogeneous expansion, which grossly fails to reproduce the evolution of coordination number and density on reducing the confining pressure (Figs. 12 and 13). The memory of larger stresses, upon decompressing, imparts wider force distributions and larger friction mobilizations in some pressure range (Fig. 4), while such reductions of coordination numbers take place. It should be expected that decompression is less predictible, because it is an evolution towards a larger disorder, and small differences can be amplified in the process. This contrasts with the compression phase, in which, for instance, the differences between configurations A and C tend to disappear. Density differences are recovered on decreasing , with the additional phenomenon that new internal states at low pressure are thus being prepared, which also differ from the initially assembled ones. While this phenomenon escapes the currently available modelling schemes, it can be noted that configurations with a high coordination number, for nearly rigid grains (low pressure or high stiffness parameter ), are extremely rare, since each contact requires a new equation to be satisfied by the set of sphere centre positions. Equilibrium states of rigid, frictionless sphere assemblies, which are the initial states for configuration series A, apart from the motion of the scarce rattlers, are isolated points in configuration space, because of isostaticity, as discussed in paper I Agnolin and Roux (2007a). As the pressure cycle, at the microscopic scale, is not reversible, due to friction and to geometric changes, one should not expect such exceptional configurations to be retrieved upon decreasing the pressure.

We thus conclude that the internal state of granular packings, in addition to the assembling process, the effect of which was studied in paper I Agnolin and Roux (2007a), varies according to the history of stress intensities, even though, unlike in cohesive materials Wolf et al. (2005); Gilabert et al. (2007), and in contrast with changes in stress directions, such loading modes only entail very small irreversible strains. Such commonly used characteristics of granular packings as coordination number, force distribution and friction mobilization level are sensitively affected by their compression history, while strains and density changes remain very small after the assembling stage. In particular, large coordination numbers associated with an ideally successful suppression of friction in the sample preparation stage seem even more unlikely to occur generally in isotropic sphere assemblies close to the RCP density, because they do not survive compression cycles. Elastic properties are studied in paper III Agnolin and Roux (2007b), where we relate them to the microstucture of such states, thereby allowing for compararisons of numerical results to experimental ones.

As possible developments of the present study, one may simulate the effects of irreversible contact deformation, due to material plasticity or particle breakage.


  1. LMSGC is a joint laboratory depending on Laboratoire Central des Ponts et Chaussées, École Nationale des Ponts et Chaussées and Centre National de la Recherche Scientifique


  1. I. Agnolin and J.-N. Roux, Internal states of model isotropic granular packings. assembling process, geometry and contact networks., first companion paper in the series (2007a).
  2. D. M. Wood, Soil Behaviour and Critical State Soil Mechanics (Cambridge University Press, 1990).
  3. J. Biarez and P.-Y. Hicher, Elementary Mechanics of Soil Behaviour (A. A. Balkema, Rotterdam, 1993).
  4. J. K. Mitchell, Fundamentals of soil behavior (Wiley, New York, 1993).
  5. H. J. Herrmann, J.-P. Hovi, and S. Luding, eds., Physics of Dry Granular Media (Balkema, Dordrecht, 1998).
  6. H. Hinrichsen and D. E. Wolf, eds., The Physics of Granular Media (Wiley-VCH, Berlin, 2004).
  7. R. García Rojo, H. J. Herrmann, and S. McNamara, eds., Powders and Grains 2005 (Balkema, Leiden, 2005).
  8. H. Makse, N. Gland, D. Johnson, and L. Schwartz, Physical Review Letters 83, 5070 (1999).
  9. C. Thornton, Géotechnique 50, 43 (2000).
  10. H. Makse, D. Johnson, and L. Schwartz, Physical Review Letters 84, 4160 (2000).
  11. L. E. Silbert, D. Ertaş, G. S. Grest, T. C. Halsey, and D. Levine, Physical Review E 65, 031304 (2002a).
  12. L. E. Silbert, G. S. Grest, and J. W. Landry, Physical Review E 66, 061303 (2002b).
  13. A. S. J. Suiker and N. A. Fleck, ASME Journal of Applied Mechanics 71, 350 (2004).
  14. A. Donev, S. Torquato, and F. H. Stillinger, PRE 71, 011105 (2005).
  15. C. O’Hern, L. E. Silbert, A. J. Liu, and S. R. Nagel, Physical Review E 68, 011306 (2003).
  16. I. Agnolin and J.-N. Roux, Internal states of model isotropic granular packings: elastic properties, second companion paper, third in the series (2007b).
  17. F. Tatsuoka, M. Sakamoto, T. Kawamura, and S. Fukushima, Soils and Foundations 26, 65 (1986).
  18. G. Reydellet and E. Clément, Physical Review Letters 86, 3308 (2001).
  19. P. A. Vermeer, in Herrmann et al. (1998), pp. 163–196.
  20. G. Gudehus, F. Darve, and I. Vardoulakis, eds., Constitutive Relations for Soils (Balkema, Rotterdam, 1984).
  21. T. G. Thomann and R. D. Hryciw, ASTM Geotechnical Testing Journal 13, 97 (1990).
  22. S. Shibuya, F. Tatsuoka, S. Teachavorasinskun, X.-J. Kong, F. Abe, Y.-S. Kim, and C.-S. Park, Soils and Foundations 32, 26 (1992).
  23. X. Jia, C. Caroli, and B. Velický, Phys. Rev. Lett. 82, 1863 (1999).
  24. X. Jia and P. Mills, in Powders and Grains 2001, edited by Y. Kishino (Swets & Zeitlinger, Lisse, 2001), pp. 105–112.
  25. R. Kuwano and R. J. Jardine, Géotechnique 52, 727 (2002).
  26. H. Geoffroy, H. di Benedetto, A. Duttine, and C. Sauzéat, in Deformation characteristics of geomaterials, edited by H. di Benedetto, T. Doanh, H. Geoffroy, and C. Sauzéat (Swets and Zeitlinger, Lisse, 2003), pp. 353–363.
  27. M. Sharifipour, C. Dano, and P.-Y. Hicher, Wave velocities in assemblies of glass beads using bender-extender elements, Proceedings of the ”Engineering Mechanics 2004” symposium of the American Society of Civil Engineering, on CD-ROM (2004).
  28. I. Agnolin, J.-N. Roux, P. Maassad, X. Jia, and P. Mills, in García Rojo et al. (2005), pp. 313–317.
  29. A. Castellanos, Advances in Physics 54, 263 (2005).
  30. K. L. Johnson, Contact Mechanics (Cambridge University Press, 1985).
  31. D. Elata and J. G. Berryman, Mechanics of Materials 24, 229 (1996).
  32. M. R. Kuhn and C. S. Chang, International Journal of Solids and Structures 43, 6026 (2006).
  33. G.d.R. MIDI, European Physical Journal E 14, 341 (2004).
  34. F. da Cruz, S. Emam, M. Prochnow, J.-N. Roux, and F. Chevoir, Physical Review E 72, 021309 (2005).
  35. O. Pouliquen, C. Cassar, Y. Forterre, P. Jop, and M. Nicolas, in García Rojo et al. (2005), pp. 859–865.
  36. P. Jop, Y. Forterre, and O. Pouliquen, Journal of Fluid Mechanics 541, 167 (2005).
  37. P. Jop, Y. Forterre, and O. Pouliquen, Nature 44, 727 (2006).
  38. H. A. Makse, N. Gland, D. L. Johnson, and L. Schwartz, Physical Review E 70, 061302 (2004).
  39. J. D. Goddard, Proc. Roy. Soc. London 430, 105 (1990).
  40. J.-N. Roux, Physical Review E 61, 6802 (2000).
  41. S. McNamara, R. García Rojo, and H. J. Herrmann, Physical Review E 72, 021304 (2005).
  42. T. C. Halsey and D. Ertaş, PRL 83, 5007 (1999).
  43. D. E. Wolf, T. Unger, D. Kadau, and L. Brendel, in García Rojo et al. (2005), pp. 525–533.
  44. F. A. Gilabert, J.-N. Roux, and A. Castellanos, Phys. Rev. E 75, 011303 (2007).
This is a comment super asjknd jkasnjk adsnkj
The feedback cannot be empty
Comments 0
The feedback cannot be empty
Add comment

You’re adding your first comment!
How to quickly get a good reply:
  • Offer a constructive comment on the author work.
  • Add helpful links to code implementation or project page.