A Contact elasticity and friction

Internal states of model isotropic granular packings. I. Assembling process, geometry and contact networks.


This is the first paper of a series of three, in which we report on numerical simulation studies of geometric and mechanical properties of static assemblies of spherical beads under an isotropic pressure. The influence of various assembling processes on packing microstructures is investigated. It is accurately checked that frictionless systems assemble in the unique random close packing (RCP) state in the low pressure limit if the compression process is fast enough, higher solid fractions corresponding to more ordered configurations with traces of crystallization. Specific properties directly related to isostaticity of the force-carrying structure in the rigid limit are discussed. With frictional grains, different preparation procedures result in quite different inner structures that cannot be classified by the sole density. If partly or completely lubricated they will assemble like frictionless ones, approaching the RCP solid fraction with a high coordination number: on the force-carrying backbone. If compressed with a realistic coefficient of friction packings stabilize in a loose state with and . And, more surprisingly, an idealized “vibration” procedure, which maintains an agitated, collisional régime up to high densities results in equally small values of while is close to the maximum value . Low coordination packings have a large proportion (10%) of rattlers – grains carrying no force – the effect of which should be accounted for on studying position correlations, and also contain a small proportion of localized “floppy modes” associated with divalent grains. Low pressure states of frictional packings retain a finite level of force indeterminacy even when assembled with the slowest compression rates simulated, except in the case when the friction coefficient tends to infinity. Different microstructures are characterized in terms of near neighbor correlations on various scales, and some comparisons with available laboratory data are reported, although values of contact coordination numbers apparently remain experimentally inaccessible.

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

I Introduction

i.1 Context and motivations

The mechanical properties of solidlike granular packings and their microscopic, grain-level origins are an active field of research in material science and condensed-matter physics Herrmann et al. (1998); Hinrichsen and Wolf (2004); García Rojo et al. (2005). Motivations are practical, originated in soil mechanics and material processing, as well as theoretical, as general approaches to the rheology of different physical systems made of particle assemblies out of thermal equilibrium Liu and Nagel (2001) are attempted.

The packing of equal-sized spherical balls is a simple model for which there is a long tradition of geometric characterization studies. Packings are usually classified according to their density or solid volume fraction , and the frequency of occurrence of some local patterns. Direct observations of packing microstructure is difficult, although it has recently benefitted from powerful imaging techniques Richard et al. (2003); Aste et al. (2004, 2005). The concept of random close packing (RCP), is often invoked Cumberland and Crawford (1987); Bideau and Dodds (1991), although some authors criticized it as ill-defined Torquato et al. (2000). It corresponds to the common observation that bead packings without any trace of crystalline order do not exceed a maximum density, , slightly below  Cumberland and Crawford (1987).

Mechanical studies in the laboratory have been performed on granular materials for decades in the realm of soil mechanics, and the importance of packing fraction on the rheological behavior has long been recognized Wood (1990); Biarez and Hicher (1993); Mitchell (1993); Herrmann et al. (1998). The anisotropy of the packing microstructure, due to the assembling process, has also been investigated Oda (1972); Arthur and Menzies (1972), and shown to influence the stress-strain behavior of test samples Fukushima and Tatsuoka (1984), as well as the stress field and the response to perturbations of gravity-stabilized sandpiles or granular layers Vanel et al. (1999).

Discrete numerical simulation Cundall and Strack (1979) proved a valuable tool to investigate the internal state of packings, as it is able to reproduce mechanical behaviors, and to identify relevant variables other than , such as coordination number and fabric (or distribution of contact orientations) Bathurst and Rothenburg (1990); Radjai and Roux (2001, 2004); Rothenburg and Kruyt (2004); Roux and Chevoir (2005). In the case of sphere packings, simulations have been used to characterize the geometry of gravity-deposited systems Silbert et al. (2002a, b) or oedometrically compressed ones Makse et al. (2000), to investigate the quasistatic, hysteretic stress-strain dependence in solid packings Thornton (2000); Suiker and Fleck (2004), and their pressure-dependent elastic moduli in a compression experiment Makse et al. (1999, 2004).

However, in spite of recent progress, quite a few basic questions remain unresolved. It is not obvious how closely the samples used in numerical simulations actually resemble laboratory ones, for which density is often the only available state parameter. Both simulations and experiments resort to certain preparation procedures to assemble granular packings, which, although their influence is recognized as important, are seldom studied, or even specified. One method to produce dense packings in simulations is to set the coefficient of intergranular friction to zero  Thornton (2000); Makse et al. (2000) or to a low value Suiker and Fleck (2004) in the assembling stage, while a granular gas gets compacted and equilibrated under pressure. On keeping frictionless contacts, this results, in the limit of low confining stress, in dense systems with rather specific properties O’Hern et al. (2003); Makse et al. (2004), related to isostaticity and potential energy minimization Roux (2000). Examples of traditional procedures in soil mechanics are rain deposition under gravity, also known as air pluviation (which produces homogeneous states if grain flow rate and height of free fall Rad and Tumay (1987) are maintained constant) and layerwise deposition and dry or moist tamping. Those two methods were observed to produce, in the case of loose sands, different structures for the same packing density Benahmed et al. (2004). Densely packed particle assemblies can also be obtained in the laboratory by vibration, or application of repeated “taps” Nowak et al. (1998); Philippe and Bideau (2002) to a loose deposit. How close are dense experimental sphere packings to model configurations obtained on simulating frictionless particles ? How do micromechanical parameters influence the packing structure ? Is the low pressure limit singular in laboratory grain packings and in what sense ?

i.2 Outline of the present study

The present paper provides some answers to such questions, from numerical simulations in the simple case of isotropically assembled and compressed homogeneous packings of spherical particles. It is the first one of a series of three, and deals with the geometric characterization of low pressure isotropic states assembled by different procedures, both without and with intergranular friction. The other two, hereafter referred to as papers II Agnolin and Roux (2007a) and III Agnolin and Roux (2007b), respectively investigate the effects of compressions and pressure cycles, and the elastic response of the different numerical packings, with comparisons to experimental results. Although mechanical aspects are hardly dealt with in the present paper, we insist that geometry and mechanics are strongly and mutually related. We focus here on the variability of the coordination numbers, which will prove important for mechanical response properties of granular packings, and show that equilibrated packs of identical beads can have a relatively large numbers of “rattler” grains, which do not participate in force transmission. We investigate the dependence of initial states on the assembling procedure, both with and without friction. We study the effects of procedures designed to produce dense states (close to RCP), and we characterize the geometry of such states on different scales.

It should be emphasized that we do not claim here to mimic experimental assembling procedures very closely. Rather, we investigate the results of several preparation methods, which are computationnally convenient, maintain isotropy, and produce equilibrated samples with rather different characteristics. Those methods nevertheless share some important features with laboratory procedures, and we shall argue that the resulting states are plausible models for experimental samples.

The numerical model and the simulation procedures (geometric and mechanical parameters, contact law, boundary conditions) are presented in Section II, where some basic definitions and mechanical properties pertaining to granular packs are also presented or recalled. Part III discusses the properties of frictionless packings, and introduces several characterization approaches used in the general case as well. Section IV then describes different assembling procedures of frictional packings and the resulting microstructures. Section V discusses perspectives to the present study, some of which are pursued in papers II and III of the series. Appendices deal with technical issues, and also present a more detailed comparison with some experimental data.

This being a long paper, it might be helpful to specify which parts can be read independently. On first going through the paper, the reader might skip Section III.4, dealing with a rather specific issue. The properties stated or recalled in Section II.3 are used to discuss stability issues and isostatic values of coordination numbers, but they can also be overlooked in a first approach. Finally, Section IV can be read independently from Section III, apart from the explanations about equilibrium conditions (in paragraph III.2.2) and the treatment of rattlers (paragraph III.5.1). Sections III and IV both have conclusive subsections which summarize the essential results.

Ii Model, numerical procedures, basic definitions

ii.1 Intergranular forces

We consider spherical beads of diameter (the value of which, as we ignore gravity, will prove irrelevant), interacting in their contacts by the Hertz law, relating the normal force to the elastic normal deflection as :


In Eqn. 1, we introduced the notation

is the Young modulus of the beads, and the Poisson ratio. For spheres, , the elastic deflection of the contact, is simply the distance of approach of the centers beyond the first contact. The normal stiffness of the contact is defined as the rate of change of the force with normal displacement:


Although many geometric features of particle packings do not depend on the details of the model for contact elasticity, and could be observed as well with a simpler, linear unilateral elastic model, it is necessary to implement suitable non-linear contact models to deal with the mechanical properties in papers II and III Agnolin and Roux (2007a, b). Tangential elasticity and friction in contacts are appropriately described by the Cattaneo-Mindlin-Deresiewicz laws Johnson (1985), which we implement in a simplified form, as used e.g., in refs. Makse et al. (1999); Somfai et al. (2005): the tangential stiffness relating, in the elastic regime, the increment of tangential reaction to the relative tangential displacement increment is a function of (or ) alone (i.e., it is kept constant, equal to its value for ):


To enforce the Coulomb condition with friction coefficient , has to be projected back onto the circle of radius in the tangential plane whenever the increment given by eqn. 3 would cause its magnitude to exceed this limit. Moreover, when decreases to , is scaled down to the value it would have had if had constantly been equal to in the past. It is not scaled up when increases. Such a procedure, suggested e.g., in Elata and Berryman (1996), avoids spurious increases of elastic energy for certain loading histories. More details are given in Appendix A.

Finally, tangential contact forces have to follow the material motion. Their magnitudes are assumed here not to be affected by rolling (i.e., rotation about a tangential axis) or pivoting (i.e., rotation about the normal axis), while their direction rotates with the normal vector due to rolling, and spins around it with the average spinning rate of the two spheres (to ensure objectivity). The corresponding equations are given in Appendix B.

In addition to the contact forces specified above, we introduce viscous ones, which oppose the normal relative displacements (we use the convention that positive normal forces are repulsive):


The damping coefficient depends on , and we choose its value as a fixed fraction of the critical damping coefficient of the normal (linear) spring of stiffness (as given by (2)) joining two beads of mass :


From (2), is thus proportional to , or to . The same damping law was used in Somfai et al. (2005). Admittedly, the dissipation given by (4)-(5) has little physical justification, and is rather motivated by computational convenience. We shall therefore assess the influence of on the numerical results. The present study being focussed on statics, we generally use a strong dissipation, , to approach equilibrium faster. This particular value is admittedly rather arbitrary: the initial motivation for choosing is the computational inefficiency of overdamped contacts with in the case of linear contact elasticity. Yet we did not check whether values of 1 or even higher would cause any problem with Hertzian contacts. In the linear case, the restitution coefficient in a binary collision varies as a very fastly decreasing function of , and changes of in the range between and 1 have virtually no detectable effect.

We do not introduce any tangential viscous force, and impose the Coulomb inequality to elastic force components only. We choose the elastic parameters  GPa and , suitable for glass beads, and the friction coefficient is attributed a moderate, plausible value . These choices are motivated by comparisons to experimental measurements of elastic moduli, to be carried out in paper III Agnolin and Roux (2007b).

ii.2 Boundary conditions and stress control

The numerical results presented below were obtained on samples of beads, enclosed in a cubic or parallelipipedic cell with periodic boundary conditions.

It is often in our opinion more convenient to use pressure (or stress) than density (or strain) as a control parameter (a point we discuss below in Section III). We therefore use a stress-controlled procedure in our simulations, which is adapted from the Parrinello-Rahman molecular dynamics (MD) scheme Parrinello and Rahman (1981). The simulation cell has a rectangular parallelipipedic shape with lengths parallel to coordinate axes (). values might vary, so that the system has configurational degrees of freedom,which are the positions and orientations of the particles and lengths . denotes the sample volume. We seek equilibrium states with set values of all three principal stresses . We use the convention that compressive stresses are positive.

It is convenient to write position vectors , defining a square matrix with ’s on the diagonal, as

denoting corresponding vectors in a cubic box of unit edge length. In addition to particle angular and linear velocities, which read

one should evaluate time derivatives . Equations of motion are written for particles in the standard form, i.e., ( denoting the total force exerted on grain )


and the usual equation for angular momentum. Meanwhile, lengths satisfy the following equation of motion, in which is the vector joining the center of to the center of , subject to the usual nearest image convention of periodic boundary conditions:


Within square brackets on the right-hand side of Eqn. 7, one recognizes the familiar formula Christoffersen et al. (1981); Hansen and McDonald (1986); Bathurst and Rothenburg (1990) for , being the average stress in the sample:


All three diagonal stress components should thus equate the prescribed values at equilibrium. The acceleration term will cause the cell to expand in the corresponding direction if the stress is too high, and to shrink if it is too low. Eqn. 7 involves a generalized mass associated with the changes of shape of the simulation cell. is set to a value of the order of the total mass of all particles in the sample. This choice was observed to result in collective degrees of freedom approaching their equilibrium values under prescribed stress somewhat more slowly (but not exceedingly so) than (rescaled) positions .

The original Parrinello-Rahman method was designed for conservative molecular systems, in such a way that the set of equations is cast in Lagrangian form. This implies in particular additional terms in (6), involving . Such terms were observed to have a negligible influence on our calculations and were consequently omitted. Granular materials are dissipative, and energy conservation is not an issue (except for some elastic properties, see paper III Agnolin and Roux (2007b)). Further discussion of the stress-controlled method is provided in Appendix C.

Equations (6) and (7), with global degrees of freedom slower than particle positions, lead to dynamics similar to those of a commonly used procedure in granular simulation Cundall (1988). This method consists in repeatedly changing the dimensions of the cell by very small amounts, then computing the motion of the grains for some interval of time. A “servo mechanism” can be used to impose stresses rather than strains Makse et al. (2004). Our approach might represent a simplification, as it avoids such a two-stage procedure. It should be kept in mind that we restricted our use of equation (7) to situations when changes in the dimensions of the simulation cell are very slow and gradual. The perturbation introduced in the motion of the grains, in comparison to the more familiar case of a fixed container, is very small.

ii.3 Rigidity and stiffness matrices

We introduce here the appropriate formalism and state the relevant properties of static contact networks. It is implied throughout this section that small displacements about an equilibrium configurations are dealt with to first order (as an infinitesimal motion, i.e. just like velocities), and related to small increments of applied forces, moments and stresses. In the following we shall exploit the definitions of stiffness matrices (Eqn. 18) and to discuss stability properties of packings. The corrections to the degree of force indeterminacy due to free mechanism motions, as expressed by relations (19) or (20), will also be used.

The properties are stated in a suitable form to the periodic boundary conditions with controlled diagonal stress components, as used in our numerical study.

Definition of stiffness matrix

We consider a given configuration with bead center positions (, ) and orientations (, ), and cell dimensions (, ). The grain center displacements are conveniently written as

with a set of displacements satisfying periodic boundary conditions in the cell with the current dimensions, and the elements of the diagonal strain matrix express the relative shrinking deformation along each direction, . Gathering all coordinates of particle (periodic) displacements and rotation increments, along with strain parameters, one defines a displacement vector in a space with dimension equal to the number of degrees of freedom ,


Let denote the number of intergranular contacts. In every contacting pair -, we arbitrarily choose a “first” grain and a “second” one . The normal unit vector points from to (along the line joining centers for spheres). The relative displacement is defined for spherical grains with radius as


in which is the vector pointing from the center of the first sphere to the nearest image of the center of the second one . The normal part of is the increment of normal deflection in the contact. (10) defines a matrix which transforms into the -dimensional vector of relative displacements at contacts :


In agreement with the literature on rigidity theory of frameworks Servatius and Servatius (1998) ( is termed normalized rigidity matrix in that reference), we call the rigidity matrix.

In each contact a force is transmitted from to , which is split into its normal and tangential components as . The static contact law (without viscous terms) expressed in Eqns. (1), (3), with the conditions stated in Section II.1, relates the -dimensional contact force increment vector , formed with the values , of the normal and tangential parts of all contact force increments, to :


This defines the matrix of contact stiffnesses . is block diagonal (it does not couple different contacts), and is conveniently written on using coordinates with as the first basis unit vector. In simple cases the block of corresponding to contact , is diagonal itself and contains stiffnesses and (twice in 3 dimensions) as given by (2) and (3):


More complicated non-diagonal forms of , which actually depend on the direction of the increments of relatives displacements in the contact, are found if friction is fully mobilized (which does not happen in well-equilibrated configurations), or corresponding to those small motions reducing the normal contact force. The effects of such terms is small, with our choice of parameters, and is discussed in paper III Agnolin and Roux (2007b).

External forces and moments (at the center) applied to the grains, and diagonal Cauchy stress components can be gathered in one -dimensional load vector :


chosen such that the work in a small motion is equal to . The equilibrium equations – the statements that contact forces balance load – is simply written with the tranposed rigidity matrix, as


This is of course easily checked on writing down all force and moment coordinates, as well as the equilibrium form of stresses:


As an example, matrices and were written down in McNamara et al. (2005) in the simple case of one mobile disk with 2 contacts with fixed objects in 2 dimensions, the authors referring to as the “contact matrix”. The same definitions and matrices are used in McNamara and Herrmann (2006) in the more general case of a packing of disks.

Returning to the case of small displacements associated with a load increment , one may write, to first order in ,


with a total stiffness matrix , comprising two parts, and , which we shall respectively refer to as the constitutive and geometric stiffness matrices. results from Eqns. 11, 12 and 15:


is due to the change of the geometry of the packing. Its elements (see Appendix B), relatively to to their counterparts in , are of order , and therefore considerably smaller in all practical cases. The constitutive stiffness matrix is also called “dynamical matrix” O’Hern et al. (2003); Somfai et al. (2005). One advantage of decomposition (18) is to separate out the effects of the contact constitutive law, contained in and those of the contact network, contained in . is sensitive in general to the orientations of normal unit vectors and to the “branch vectors” joining the grain centers to contact positions – which reduce to for spheres of radius . , on the other hand, unlike , is sensitive to the curvature of grain surfaces at the contact point Kuhn and Chang (2006); Bagi (2007).

Properties of the rigidity matrix

To the rigidity matrix are associated the concepts (familiar in structural mechanics) of force and velocity (or displacement) indeterminacy, of relative displacement compatibility and of static admissibility of contact forces. Definitions and properties stated in Roux (2000) for frictionless grains, straightforwardly generalize to packings with friction.

The degree of displacement indeterminacy (also called degree of hypostaticity Roux (2000)) is the dimension of the kernel of , the elements of which are displacements vectors which do not create relative displacements in the contacts: . Such displacements are termed (first-order) mechanisms. Depending on boundary conditions, a grain packing might have a small number of “trivial” mechanisms, for which the whole system moves as one rigid body. In our case, attributing common values of to all grains gives independent global rigid motions.

The degree of force indeterminacy h (also called degree of hyperstaticity Roux (2000)) is the dimension of the kernel of , or the number of independent self-balanced contact force vectors. If the coordinates of are regarded as the unknowns in system of equations (15), and if is supportable, then there exists a whole h-dimensional affine space of solutions.

From elementary theorems in linear algebra one deduces a general relation between h and  Roux (2000)


An isostatic packing is defined as one devoid of force and velocity indeterminacy (apart from trivial mechanisms). Excluding trivial mechanisms (thus reducing to ), and loads that are not orthogonal to them, one then has a square, invertible rigidity matrix. To any load corresponds a unique set of equilibrium contact forces. To any vector of relative contact displacements corresponds a unique displacement vector.

With frictionless objects, in which contacts only carry normal forces, it is appropriate to use -dimensional contact force and relative displacement vectors, containing only normal components, and to define the rigidity matrix accordingly Roux (2000). Then (19) should be written as


In the case of frictionless spherical particles, all rotations are mechanisms, hence a contribution of 3 to . Thus one may in addition ignore all rotations, and subtract both from and from , so that (20) is still valid. In such a case, the rigidity matrix coincides (up to a sign convention and normalization of its elements) with the one introduced in central-force networks, trusses and tensegrity structures Thorpe and Duxbury (1998). Donev et al., in a recent publication on sphere packings Donev et al. (2005), call rigidity matrix what we defined as its transpose .

ii.4 Control parameters

The geometry and the mechanical properties of sphere packings under given pressure depends on a small set of control parameters, which can be conveniently defined in dimensionless form Combe and Roux (2003); Roux and Chevoir (2005).

Such parameters include friction coefficient and viscous dissipation parameter , which were introduced in Sec. II.1.

The elastic contact law introduces a dimensionless stiffness parameter , which we define as:


Note that does not depend on bead diameter . Under pressure , the typical force in a contact is of order . It corresponds to a normal deflection such that due to the Hertz law (1). Therefore, sets the scale of the typical normal deflection in Hertzian contacts, as .

In the case of monodisperse sphere packings in equilibrium in uniform state of stress , pressure is directly related to the average normal force . Let us denote as the solid fraction and the coordination number (). As a simple consequence of the classical formula for stresses recalled in Sec. II.2 (Eqn. 8 in the static case, or Eqn. 16), one has, neglecting contact deflections before diameter ,


whence an exact relation between and contact deflections:

The limit of rigid grains is approached as . can reach very high values for samples under their own weight, but most laboratory results correspond to levels of confining pressure in the 100 kPa range. Experimental data on the mechanical properties of granular materials in quasistatic conditions below a few tens of kPa are very scarce (see, however, Tatsuoka et al. (1986) and Reydellet and Clément (2001)). This is motivated by engineering applications (100 kPa is the pressure below a few meters of earth), and this also results from difficulties with low confining stresses. Below this pressure range, stress fields are no longer uniform, due to the influence of the sample weight, and measurements are difficult (e.g., elastic waves of measurable amplitude are very strongly damped).

We set the lowest pressure level for our simulation of glass beads to 1 kPa or 10 kPa, which corresponds to and . Such values, as we shall check, are high enough for some characteristic properties of rigid sphere packings to be approached with good accuracy. Upon increasing , the entire experimental pressure range will be explored in the two companion papers Agnolin and Roux (2007a, b).

Another parameter associated with contact elasticity is the ratio of tangential to normal stiffnesses (constant in our model), related to the Poisson ratio of the material the grains are made of. Although we did not investigate the role of this parameter, several numerical studies Combe (2002); Somfai et al. (2005) showed its influence on global properties to be very small.

The “mass” of the global degrees of freedom is chosen to ensure slow and gradual changes in cell dimensions, and dynamical effects are consequently 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 . was successfully used as a control parameter in dense granular shear flows MIDI (2004); da Cruz et al. (2005); Pouliquen et al. (2005), which might be modelled on writing down the dependence of internal friction and density Jop et al. (2005, 2006).

The sensitivity to dynamical parameters and should be larger in the assembling stage (as studied in the present paper) than in the subsequent isotropic compression of solid samples studied in paper II Agnolin and Roux (2007a), for which one attempts to approach the quasi-static limit. In the following we will assess the influence of parameters , , and on sample states and properties.

Iii Low-pressure isotropic states of frictionless packings.

iii.1 Motivations

Numerical samples are most often produced by compression of an initially loose configuration (a granular gas) in which the grains do not touch. If the friction coefficient is set to zero at this stage, one obtains dense samples, which depend very little on chosen mechanical parameters. These frictionless configurations are in a particular reference state which was recently investigated by several groups O’Hern et al. (2003); Donev et al. (2005). We shall dwell on such an academic model as assemblies of rigid or slightly deformable frictionless spheres in mechanical equilibrium for several reasons. First, we have to introduce various characterizations of the microstructure of sphere packings that will be useful in the presence of friction too. Then, such systems possess rather specific properties, which are worth recalling in order to assess whether some of them could be of relevance in the general case. Frictionless packings also represent, as we shall explain, an interesting limit case. Finally, one of our objectives is to establish the basic uniqueness, in the statistical sense, of the internal state of such packings under isotropic, uniform pressures, provided crystallization is thwarted by a fast enough dissipation of kinetic energy.

iii.2 Assembling procedures

Previous results

Since we wish to discuss a uniqueness property, we shall compare our results to published ones whenever they are available. Specifically, we shall repeatedly refer to the works of O’Hern, Silbert, Liu and Nagel O’Hern et al. (2003), and of Donev, Torquato and Stillinger Donev et al. (2005), hereafter respectively abbreviated as OSLN and DTS. Both are numerical studies of frictionless sphere packings under isotropic pressures.

OSLN use elastic spheres, with either Hertzian or linear contact elasticity. They control the solid fraction , and record the pressure at equilibrium. Their samples (from a few tens to about 1000 spheres) are requested to minimize elastic energy at constant density. For each one, pressure and elastic energy vanish below a certain threshold packing fraction , which is identified to the classical random close packing density. Above , pressure and elastic constants are growing functions of density. OSLN report several power law dependences of geometric and mechanical properties on which we shall partly review.

DTS differ in their approach, as unlike OSLN (and unlike us) they use strictly rigid spherical balls, and approach the density of equilibrated rigid, frictionless sphere packing from below. They use a variant of the classical (event-driven) hard-sphere molecular dynamics method Alder and Wainwright (1959); Hansen and McDonald (1986), in which sphere diameters are continuously growing, the Lubachevsky-Stillinger (LS) algorithm Lubachevsky and Stillinger (1990); Lubachevsky et al. (1991), to compress the samples. DTS’s approximation of the strictly rigid sphere packing as the limit of a hard sphere glass with very narrow interstices (gaps) between colliding neighbors (contact forces in the static packing are then replaced by transfers of momentum between neighbors), and their resorting to linear optimization methods Roux (2000); Donev et al. (2004), enable them to obtain very accurate results in samples of 1000 and 10000 beads.

DTS expressed doubts as to whether numerical “soft” (elastic) sphere systems could approach the ideal rigid packing properties, and both groups differ in their actual definition of jamming and on the relevance and definition of the random close packing concept. Relying on our own simulation results, we shall briefly discuss those issues in the following.

Frictionless samples obtained by MD

Our numerical results on packings assembled without friction are based on five different configurations of beads prepared by compression of a granular gas without friction. First, spheres are placed on the sites of an FCC lattice at packing fraction (below the freezing density,  Volkov et al. (2002)). Then they are set in motion with random velocities, and left to interact in collisions that preserve kinetic energy, just like the molecules of the hard-sphere model fluid studied in liquid state theory Hansen and McDonald (1986); Volkov et al. (2002). We use the traditional event-driven method Alder and Wainwright (1959), in a cubic cell of fixed size, until the initial crystalline arrangement has melted. Then, velocities are set to zero, and the molecular dynamics method of Section II is implemented with an external pressure equal to 10 kPa for glass beads (). Energy is dissipated thanks to viscous forces in contacts, and the packing approaches an equilibrium state. Calculations are stopped when the net elastic force on each particle is below , the elastic contributions to the stress components equal the prescribed value with relative error smaller than and the kinetic energy per particle is below . On setting all velocities to zero, it is observed that the sample does not regain kinetic energy beyond that value, while the unbalanced force level does not increase. We have thus a stable equilibrium state. This is further confirmed by the absence of mechanism in the force-carrying contact network, apart from the trivial free translational motion of the whole set of grains as one rigid body. From (18), mechanisms coincide with “floppy modes” of the constitutive stiffness matrix (i.e., the elements of its kernel). The geometric stiffness , as checked in Appendix B, is a very small correction (compared to those of the elements of matrix are of order ).

In the following, such configurations assembled without friction will be referred to as A states.

In order to check for a possible influence of the assembling procedure on the final configurations, we simulated another, similar sample series, denoted as A’, for which the LS algorithm was used to bring the solid fraction from to , before equilibrating at the desired pressure with Hertzian sphere molecular dynamics.

Observed geometric and mechanical characteristics of A and A’ states are reported below and compared to other published results, in particuler those of OSLN and DTS. We also state specific properties of rigid frictionless sphere packings, to which A configurations at high are close. Unlike OSLN, we use pressure or stiffness level as the control parameter. The state OSLN refer to as “point J”, which appears as a rigidity threshold if solid fraction is used as the control parameter, is approached here as .

Compression rates and duration of agitation stage

Molecular dynamics is not the fastest conceivable route to minimize the sum of elastic and potential energies, and the MD approach does not necessarily find the nearest minimum in configuration space. For that purpose, the direct conjugate gradient minimization approach, as used by OSLN, which involves no inertia and follows a path of strictly decreasing energy in configuration space is the best candidate.

However, the time scales involved in the MD simulations can be compared to experimental ones. In simulations, A configurations approach their final density within a few tens of time , and come to their final equilibrium with a few hundreds of . Comparable laboratory experiments in which dense samples are assembled are sample preparation with the pluviation or rain deposition technique, in which grains are deposited at constant flow rate under gravity, with a constant height of free fall Rad and Tumay (1987); Emam et al. (2005); Benahmed (2001); Benahmed et al. (2004). Such an assembling technique produces homogeneous samples. Grains are first agitated near the free surface, and then subjected to a quasistatic pressure increase as pouring procedes. The relevant pressure scale corresponds therefore to the weight of the agitated superficial layer of the sample being assembled Emam et al. (2005), typically of the order of 10 diameters, hence and , about  s for  mm. Approximating the compaction time by the time needed to renew entirely the agitated superficial layer, we obtain a few times  s if this time is to be of order , as in our simulations. This corresponds to a fraction of a second to fill up a  cm high container, a value within the experimental range. The main conclusion from this crude analysis is that laboratory assembling processes are rather fast, with typical compaction times similar to those of our simulated isotropic compression procedure.

On the other hand, the LS procedure followed by DTS, which we used to produce our A’ samples, unavoidably involves many collisions and a considerable level of agitation while particle diameters grow at a prescribed rate. In practice, kinetic energy actually increases on implementing the LS algorithm: receding velocities after a collision have to be artificially increased in order to make sure particles that are continuously growing in diameter actually move apart after colliding Lubachevsky and Stillinger (1990); Lubachevsky et al. (1991). Velocities have to be scaled down now and then for computational convenience, a feature the actual compacting process, depending on the ratio of growing rate to quadratic velocity average, is sensitive to. In our implementation there were typically 110 collisions per sphere in the range (the most dangerous interval for crystal nucleation Volkov et al. (2002); Luchnikov et al. (2002)), and 90 collisions for . DTS report using expansion rates of , while ours started as , in units of the quadratic mean velocity.

Consequently, the order of simulation results, from the fastest, least agitated case to the slowest one is as follows: first the OSLN results, then our A, followed by our A’ series, and finally the simulations by DTS (who used a slower implementation of the LS method than our A’ one).

iii.3 Energy minimization and density

What is “jamming” ?

In spite of a long tradition of studies on the geometry of sphere assemblies, the connection between mechanical equilibrium and density maximization has seldom been stressed. This property was presented, in slightly different forms, in the mathematics Connelly (1988) and physics Roux (2000); Donev et al. (2005) literature. It is worth recalling it here, as the purpose of this work is to discuss both geometric and mechanical properties of such particle packings. This connection is simply expressed on noting that configurations of rigid, frictionless, non-adhesive spherical particles in stable equilibrium under an isotropic confining pressure are those that realize a local minimum of volume in configuration space, under the constraint of mutual impenetrability. It is no wonder then that the isotropic compaction of frictionless balls is often used as a route to obtain dense samples Thornton (2000); Makse et al. (2004). In DTS Donev et al. (2005) and in other works by the same group Kansal et al. (2002); Donev et al. (2004), the authors use a definition of strictly jammed configurations of hard particles as those for which particles cannot move without interpenetrating or increasing the volume of the whole system. Their definition is therefore exactly equivalent to that of a stable equilibrium state with rigid, frictionless grains under an isotropic confining pressure.

If we now turn to elastic, rather than rigid, spherical particles, with Hertzian contacts as defined in Sec. II, then stable equilibrium states under given pressure are local minima of the potential energy defined as ( denotes the Heavside step function)


As stiffness parameter increases, the second term of (24) imparts an increasing energetic cost to elastic deflections , and the solution becomes an approximation to a minimum of the first term, with impenetrability constraints, i.e., a stable equilibrium state of rigid, frictionless balls. The value of is an indicator of the distance to the ideal, rigid particle configuration, and it is arguably more convenient to use that the density, used by OSLN, because it does not vary between different samples. OSLN had to adjust the density separately for each sample in order to approach the limit of rigid grains, so that the pressure approached zero, corresponding to a rigidity threshold. Their definition of jamming is based on a local minimum of elastic energy, and therefore also coincides with ours: a jammed state is a stable equilibrium state.

Solid fractions

Our A configurations have a solid fraction (indicated error bars correspond throughout the paper to one sample-to-sample standard deviation). We shall check below that the small density difference between and is much smaller than the statistical uncertainty on . OSLN performed a careful statistical analysis of finite size effects and uncertainty on , leading to estimates shown on Fig. 1. Fig. 1 also shows another MD data point we obtained for . Our values coincide with OSLN’s estimation of size-dependent averages and fluctuations, once it is extrapolated to larger sample sizes (or, possibly, our configurations are very slightly denser). Our A’ samples exhibit higher densities than A ones, – a fairly small difference, but clearly larger than error bars.

Figure 1: (Color online) Solid fraction versus . Blue dotted line: average value, with standard deviation indicated as error bars, according to OSLN’s results, extrapolated to . Black round dots with error bars: our A samples for and other, very similar results for . Square dot: A’ samples with (this point has a smaller error bar).

DTS do not report values very precisely, but mention solid fractions in the range 0.625 to 0.63 (Donev et al., 2005, page 7), on excluding the volume of rattlers, particles that transmit no force. This entails once those inactive grains, which represent about 2.2% of the total number, are taken into account. LS-made samples were shown in Torquato et al. (2000) to jam, depending on the compression rate, over the whole solid fraction range between 0.64 and the maximum value corresponding to the perfect FCC crystal. The final values of the solid fraction therefore correlate with the duration of the agitation stage in the assembling procedure. The RCP density is traditionnally associated with a minimization of crystalline order. In the next section we check for indications of incipient crystalline order in A and A’ samples.

iii.4 Traces of crystalline order

The possible presence of crystal nuclei, the FCC and HCP lattices (the former the more stable thermodynamically) and hybrids thereof being the densest possible arrangements, is a recurring issue in sphere packing studies. A recent numerical study of crystallization dynamics in the hard sphere fluid is that of Volkov et al. Volkov et al. (2002), in which the authors used several indicators and measures of incipient crystalline order that we apply here to A and A’ states. First, bonds are defined as (fictitious) junctions between the centers of neighboring spheres if their distance is smaller than some threshold, often chosen as corresponding to the first minimum in pair corelation function (about in our case, see Sec. III.6.1). Then, a local order parameter is associated to each grain , as:


in which is an average over all neighbors of numbered from to , the number of bonds of :


denoting as usual the unit vector pointing from the center of to the center of .

and , in particular, have been used to distinguish different local orders Volkov et al. (2002); Aste et al. (2005); Luchnikov et al. (2002). In the sequel we use the average of (26) over all grains, as well as a global parameter , defined on taking the average over all bonds within the sample, instead of those of a particular grain in (25). The values of those parameters are given in table 1. Global values are small in large samples, because they tend to average to zero in the presence of randomly oriented polycrystalline textures. They can be used nevertheless to observe crystallization in samples of beads, as they finally reach values comparable to the perfect crystal one Luchnikov et al. (2002).

Next, following Volkov et al. (2002), we normalize the set , on multiplying, for any given and , each of its components by an appropriate common factor thus obtaining , such that


If the values are viewed as the components of a -dimensional local order parameter, then might be viewed as its “phase” or “angular” part, characteristic of the choice of a direction, rather than of the intensity or extent with which the system is locally ordered. Then a bond is termed crystalline if it joins two particles for which those “phases” are sufficiently correlated: (the star indicates complex conjugation)


A particle is said to be in a crystalline configuration if at least 7 of its bonds (out of 12.5-13, see table 1)) is “crystalline”, according to definition (28) with . One may check how numerous those particles are and whether they tend to cluster in crystalline regions. Table 1 contains those various indicators, as observed in samples of type A and A’ at the largest studied stiffness level. Order parameters have a very small value, indicating as expected a large distance to crystal order. Only a small fraction of bonds and grains are declared “crystalline” according to the above definitions. However, it does transpire from the data of table 1 that A’ states are consistently more “ordered” than A ones, with a small, but systematic difference for all listed indicators (see also Appendix D).

Table 1: Indicators of possible incipient crystalline order in states A and A’ at . is the coordination number of first heighbors, the “crystalline bond” coordination number, and the global and (average) local order parameters, its average value within crystalline regions, the fraction of “crystalline” particles and the mass average of the number of particles in a “crystalline cluster”. First neighbors are defined here as those closer than distance or , near the first minimum in .

Most notable is the increase of the size of “crystalline” regions. A direct visualization of those domains, as we checked, shows that they are quite far from perfectly ordered, but reveals some local tendency to organization in parallel stacked layers, and to the formation of 2D triangular lattice patterns within the layers. Luchnikov et al. Luchnikov et al. (2002) report simulation of 16000 particle samples of the hard sphere fluid evolving towards crystallization at constant density (between and ), as monitored by the global parameter and the distribution of local values. They observed, like Volkov et al. Volkov et al. (2002), that several thousands of collisions per particle were necessary for a significant evolution to take place, which is compatible with our observation of a detectable, but very small tendency with about 100 collisions per particle with our A’ samples.

and , as defined in (25), were also used by Aste et al. Aste et al. (2005) to characterize local arrangements, in an experimental study of sphere packing geometry by X-ray tomography. These results rely on observations of large samples of tens to hundreds of thousands of beads, although not isotropic. Particles are classified according to the pair of values , . We compared the geometry of our numerical samples of similar density to those experimental data, with the result that although the most frequently observed values of and were quite close to experimental ones in dense samples, and the proportion of hcp-like particles were similar, fcc-like local environments were exceptional in simulations, whereas a few percent of the spheres were classified in that category in the experimental results. Quantitative results are given in Appendix D. Nucleation of crystalline order is strongly sensitive to sample history and boundary conditions Volkov et al. (2002).

iii.5 Properties of force networks

Identification and treatment of rattlers

The rattlers are defined as the grains that do not participate in carrying forces and remain, therefore, free to “rattle” within the cage formed by their force-carrying, rigidly fixed neighbors. We refer to the network of contacting grains that carry forces as the backbone. The backbone is the structure formed by non-rattler grains. The fraction of rattlers at is in A samples, and it is slightly higher, in A’ ones. DTS report , and hence once again our A’ results are closer to theirs. The proportion of rattlers increases slightly for stiffer contacts (higher values).

Distinguishing between the backbone and the rattlers requires some care, as very small forces on the backbone might be confused with forces below tolerance between rattler and backbone grains, or between two rattlers. We apply the following simple procedure. First, we regard as rattlers all spheres having less than four contacts: less than three contacts implies a mechanism, and only three is impossible if forces are all strictly compressive. We also discard from the backbone all spheres with only forces smaller than the tolerance. Then, all the contacts of eliminated spheres being also removed, other spheres might (although this is an extremely rare occurrence) have less than four contacts, so the procedure is iterated (twice at most is enough in our samples, although one such sweep is usually enough) until no more rattler is detected. We found this method to work correctly for and . If one eliminates too many particles, the identified backbone might become floppy (hypostatic). We check, however, that its constitutive stiffness matrix remains positive definite, thereby avoiding such pitfall. The proportion of rattlers is likely to increase for stiffer contacts (higher ).

The presence of rattlers complicates the analysis of geometric properties of static packings, because their positions are not determined by the equilibrium requirement. The rattlers are free to move within a “cage” formed by their backbone neighbors, and there is no obvious way in principle to prefer one or another of their infinitely many possible positions. This renders the evaluation of geometric data like pair correlations somewhat ambiguous. Moreover, rattlers, although scarce in frictionless packings, can be considerably more numerous in frictional ones (see Section IV). We therefore specify whether the results correspond to direct measurements on the configurations resulting from the simulations, with rattlers floating in some positions resulting from compaction dynamics, or whether rattlers have been fixed, each one having three contacts with the backbone (or some previously fixed other rattler). To compute such fixed rattler positions with MD, we regard each backbone grain as a fixed object, exert small isotropically distributed random forces on all rattlers and let them move to a final equilibrium position (assuming frictionless contacts). A third possibility is to eliminate rattlers altogether before recording geometric data. These are three choices referred to as I, II and III in the sequel, and we denote observed quantities with superscripts I, II or III accordingly.

Packings under gravity, if locally in an isotropic state of stress, are expected to be in the same internal state and to exhibit the same properties as the ones that are simulated here. In such a situation, individual grain weights are locally, within an approximately homogeneous subsystem, dominated by the externally imposed isotropic pressure. There is no rattler under gravity, but some grains are simply feeling their own weight, or perhaps that of one or a few other grains relying on them. Such grains are those that would be rattlers in the absence of gravity. Instead of freely floating within the cage of their backbone neighbors, they are supported by the cage floor. The situation should therefore be similar to that of our samples after all rattlers have been put in contact with the force-carrying structure (treatment II), except that the small external forces applied to them are all directed downwards.

Coordination numbers

Table 2 gives the distribution of local coordination number values among the spheres for A and A’ states at . In this table, is the proportion of grains with contacts.

A (I)
A’ (I)
A (II)
A’ (II)
Table 2: Percentages of grains having contacts in A and A’ configurations, on ignoring the rare contacts with or between rattlers (I), or on fixing the rattlers onto the backbone with small (randomly oriented) forces (II).

If rattlers are stuck to the backbone (method II), one records slightly changed proportions of spheres with contacts, to which values observed within samples under gravity should be compared. Distributions of local coordination numbers observed by DTS coincide to the data of Table 2 within 1%. We attribute this small difference to the influence of contact deflections of order in the MD results, while the DTS results are closer to ideally rigid packings (approached as open gaps tend to zero).


We now discuss how the isostaticity property Ouaguenouni and Roux (1997a); Roux (1997); Moukarzel (1998a, b); Tkachenko and Witten (1999); Roux (2000); Moukarzel (2001, 2004) of equilibrium states of rigid, frictionless spheres, influences high configurations of type A.

Isostaticity is a property of the backbone, i.e. the force-carrying contact network, in equilibrium packings of rigid, frictionless spheres. It means that such networks are both devoid of hyperstaticity (force indeterminacy) and of hypostaticity (displacement indeterminacy), apart from possible trivial motions in which all force-carrying grains move as one rigid body. These two properties have different origins Roux (2000), and are not valid under the same assumptions. The absence of hyperstaticity ( with the notations of Sec. II.3) results from the generic disorder of the packing geometry. It would hold true for arbitrarily shaped rigid particles interacting by purely normal contact forces whatever the sign of those forces, and it applies to the whole packing, whatever the contacts the rattlers might accidentally have with the backbone. The absence of hypostaticity property (except for trivial mechanisms, ), on the other hand, is only guaranteed for spherical particles with compressive forces in the contacts, and it applies to the sole backbone.

Due to the isostaticity property, the coordination number should be equal to 6 in the rigid limit on the backbone. If is the number of force-carrying contacts, then the global (mechanical) coordination number is (possible contacts of the rattlers are discarded), and the backbone coordination number is defined as . , rather than , has the limit as . In A samples () one has (and hence ), the excess over the limit resulting from contacts that should open on further decreasing the pressure.

The isostaticity property can be used to evaluate the density increase due to finite particle stiffness. To first order in the small displacements between (or ) and the current finite pressure state A, one might use the theorem of virtual work Roux (2000), with the displacements that bring all overlaps to zero, and the current contact forces. Such motions leading to a simultaneous opening () of all contacts are only possible on networks with no hyperstaticity, because there is no compatibility condition on relative normal displacements Roux (2000). This yields an estimate of the increase of the solid fraction over its value in the rigid limit, as


This equality can be rearranged using the Hertz contact law (1) to relate to , and relation (22). We denote as the moment of order of the distribution of normal forces , normalized by the average over all contacts:


(29) can be rewritten as:


In the isostatic limit which is approached at large , the force distribution and its moments are determined by the network geometry, and we observed . Taking for and the values at the highest studied stiffness level (), this enables us to evaluate the density change between those configurations and the rigid limit as . As announced before this is smaller than the statistical uncertainty on , and hence this does not improve our estimation of the solid fraction of the packing of rigid particles (). Recalling , (31) means that the macroscopic relation between pressure and density has the same power law form () as the contact law (). This was observed by OSLN. It would hold, because of the isostaticity property in the rigid limit, for whatever exponent in the contact law, the prefactor of the macroscopic relation involving , a moment of the geometrically determined force distribution.

As a consequence of (31), one can simply relate the bulk modulus of frictionless packings to the pressure, as observed by OSLN too, a property which will be used and discussed in paper III Agnolin and Roux (2007b), which deals with elasticity of packings.

Force distribution

The force distribution we observe in A samples at high values approaches the one of a rigid packing, which due to isostaticity is a purely geometric property. It is represented on Figure 2.

Figure 2: (Color online) Probability distribution function of normalized contact forces in A and A’ configurations at high . The dashed line is the fit proposed by DTS: .

The data presented here are averaged over 5 samples. Because of relation (22), all samples prepared at the same pressure have the same average force, and this restores the “self-averaging” property, which OSLN observed to be lacking on using solid fraction instead of pressure as the control parameter. The choice of as a state variable, because of the finite size of the sample causing fluctuations of the threshold where vanishes, is less convenient in that respect.

Fig. 2 also shows that the form proposed by DTS to fit their data is in very good agreement with our results, except perhaps for large forces, for which it is a better fit for A’ data – thus providing additional evidence that A’ samples are closer than A ones to the DTS results.

iii.6 Geometric characterization

Pair correlation function

Pair correlations should preferably be measured either with method I or method II, as there is no reason to eliminate rattlers before studying geometric properties. Comparisons between pair correlation functions and (Fig. 3) show very little difference on the scale of one particle diameter. Results of Fig. 3 are very similar to other published ones (e.g., in DTS), with an apparent divergence as and a split second peak, with sharp maxima at and .

Figure 3: Pair correlation functions and versus in A samples at . Both definitions coincide on this scale (only the peak for is slightly different).

The pair correlation function should contain a Dirac mass at in the limit of rigid grains, which broadens into a sharp peak for finite contact stiffness. The weight of this Dirac term or peak in the neighbor intercenter distance probability distribution function is coordination number , and the shape of the left shoulder of the peak at finite is directly related to the force distribution :

This explains the observation by OSLN O’Hern et al. (2003) of the width of the peak decreasing approximatively as , while its height increases as , as the threshold density is approached from above. The form of the distribution of contact forces, which is determined by the geometry of the isostatic backbone, remains exactly the same for all small enough values of , with a scale factor proportional to , due to Eqns. 22 and 31.

The sharp drop of at and was found by DTS to go to a discontinuity in the rigid particle limit. This can be understood as follows. Each sphere has a number of first contact neighbors ( on average) at distance if the grains are rigid, and a number of second contact neighbors (i.e. particles not in contact with it, but having a contact with at least one of its first contact neighbors). Such second contact neighbors will make up for a significant fraction of particles with their centers at a distance , but none of them can be farther away. Futhermore, this leads to a systematic depletion of the corona (with ) by steric exclusion.

Near neighbor correlations

As , pair correlations are conveniently expressed with the gap-dependent coordination number . is the average number of neighbors of one sphere separated by an interstice narrower than . is the usual contact coordination number . Function has three possible different definitions , and according to the treatment of rattlers. All three of them were observed to grow as for smaller than about , constant taking slightly different values for , and . is equal to for , and is very well fitted with the value found by DTS (Donev et al., 2005, Fig. 8). deviates from this power law dependence corresponding to the rigid limit for small , of the order of the typical overlap , as shown on Fig. 4. This power law corresponds to diverging as as .

Figure 4: (Color online) Coordination number of near neighbors function of interstice . The power law regime extends to smaller and smaller values as increases.

Silbert et al., in a recent numerical study of states with high levels of rigidity Silbert et al. (2006) (), report observing to grow with an exponent closer to , although somewhat dependent on the choice of the interval for the fit. However this does not contradict our main conclusion that different numerical approaches track the same RCP state in the rigid limit.

Other properties of contact networks

Ignoring rattlers (method III) one may record the density of specific particle arrangements in the backbone. We thus find the contact network (joining all centers of interacting particles by an edge) to comprise a number of equilateral triangles, such that on average each backbone grain belongs to triangles. In the rigid limit this gives a Dirac term for in the distribution of angles between pairs of contacts of the same grain. Tetrahedra are however very scarce (as observed by DTS), involving about of the beads, and pairs of tetrahedra with a common triangle are exceptional (5 such pairs in 5 samples of 4000 beads). Pairs of triangles sharing a common base are present with a finite density, which explains the discontinuous drop at of , this being the largest possible distance for such a population of neighbor pairs.

iii.7 Conclusions

We summarize here the essential results of Section III, about frictionless packings.

Uniqueness of the RCP state

Our numerical evidence makes a strong case in favor of the uniqueness of the simulated rigid packing state made with frictionless spheres under isotropic pressure. Specifically, we observed quantitative agreements with other published results O’Hern et al. (2003); Donev et al. (2005) in the coordination numbers, the force distributions, the pair correlations and the frequency of occurrence of local contact patterns, even though different numerical methods have been used. The small remaining differences in solid fraction, proportion of rattlers, and probability of large contact forces all correlate with the duration of agitated assembling stage, which can be measured in terms of numbers of collisions per grain at a given density. This duration directly correlates to the packing fraction and to the small amount of crystalline order in the samples. We therefore checked in an accurate, quantitative way the traditional views about random close packing (RCP). The RCP state can be defined in practice as the unique state in which rigid frictionless spherical beads assemble in a static equilibrium state under isotropic pressure, in the limit of fast compression, so that the slow evolution towards crystallization has a negligible influence. The Lubachevsky-Stillinger algorithm tends to produce packings with a small but notable crystalline fraction.

Relevance of MD simulations, role of micromechanical parameters

The uniqueness of the RCP implies that dynamical parameters and have no influence on the frictionless packing structures, at least in the limit of fast compression rates.

Standard MD methods compare well with specifically designed methods that deal with rigid particles, and prove able to approach the rigid limit with satisfactory, if admittedly smaller, accuracy. Recalling that corresponds to glass beads under  kPa, it seems that laboratory samples under usual conditions might in principle (if friction mobilization can be suppressed) approach the ideal (rigid particle) RCP state.

Moreover, the time scales to assemble samples in MD simulations, if compared to estimated preparation times in the laboratory with such techniques as controlled pluviation, has the right order of magnitude. This means that the assembling proces is rather fast in experimental practice when grains are deposited under gravity, which explains why densities above RCP are not directly obtained. Of course, in practice, many procedures produce anisotropic states. Anisotropic packings of rigid, frictionless balls, under other confining stresses than a hydrostatic pressure, should differ from the RCP state, and the numerical simulations of gravity deposited packings of frictionless beads of refs. Silbert et al. (2002a, b) could be analysed in this respect. We chose here to study ideal preparation methods, and we only deal with isotropic systems.

Approach to isostaticity in the rigid limit

We checked that bead packings under typical laboratory pressures such as 10 kPa might closely approach the isostaticity property of rigid frictionless packings. We showed that some observations made by O’Hern et al. O’Hern et al. (2003) on pressure or bulk modulus dependence on density, and on the shape of the first peak of were direct consequences of this remarkable property.

Iv Low-pressure states of frictional packings obtained by different procedures

iv.1 Introduction

It is well known that the introduction of friction in granular packings tends in general to reduce density and coordination number, as observed in many recent numerical simulations (see e.g. Thornton (2000); Makse et al. (2000); Silbert et al. (2002a, b)), and that frictional granular assemblies, unlike frictionless ones, can be prepared in quite a large variety of different states. In the field of soil mechanics, sand samples are traditionnally classified by their density Wood (1990); Biarez and Hicher (1993); Mitchell (1993), which determines behaviors that have been observed in simulations of model systems as well Thornton (2000); Radjai and Roux (2004). (Inherent anisotropy of the fabric, i.,e., the one due to the assembling process, rather than induced by anisotropic stresses, is a secondary, less influential state variable Oda (1972); Arthur and Menzies (1972)). Engineering studies on sands usually resort to a conventional definition of minimum and maximum densities, based on standardized procedures emi (1992).

The motivation of the present study is to explore the range of accessible packing states, as obtained by different numerical procedures that produce homogeneous and isotropic periodic samples. We therefore chose to bypass the painstaking computations needed to mimic actual laboratory assembling methods, but we argue that our procedures produce plausible structures with similar properties.

One key result is that density alone does not determine the internal state of an isotropic packing, because the coordination number can vary independently.

The A-type configurations obtained without friction in Section IV.2 are local density maxima in configuration space (see Sec. III.3.1). Hence compaction methods can be regarded as strategies to circumvent the mobilization of intergranular friction forces. Two such procedures are studied here, in a simplified, idealized form: lubrication and vibration. We also simulate, as a reference, a state which can be regarded as a loose packing limit, at least with a definition relative to one assembling method and friction coefficient ; and we prepare, as an interesting limit from a theoretical point of view, infinite friction samples.

Assembling procedures are described in Section IV.2, geometric aspects are studied in Section IV.3, and contact network properties in Section IV.4. Section IV.5 summarizes the results.

iv.2 Assembling processes for frictional grains

Just like in the frictionless case, for each one of the packing states, we prepare 5 samples of 4000 beads, over which results are averaged, error bars corresponding to sample-to-sample fluctuations. The equilibrium criteria are those of Section III.2.2, supplemented with a similar condition on moments. To identify rattlers, we use the procedure defined in Section III.5.1, which is adapted to the case of frictional grains: spheres with as few as two contacts may carry forces (even large ones, as we shall see) and should not be regarded as rattlers.

Looser packings compressed with final friction coefficient

We used direct compression of a granular gas in the presence of friction (), another standard numerical procedure Makse et al. (2000); Thornton (2000); Suiker and Fleck (2004), in which the obtained density and coordination number are decreasing functions of  Thornton (2000); Makse et al. (2000); Silbert et al. (2002a, b). This produces rather loose samples hereafter referred to as D (B and C ones, to be defined further, denoting denser ones, closer to A, but arguably more “realistic”). D samples were made with exactly the same method as A ones (see Sec. III.2.2), except that the friction coefficient was used instead of . In principle, D configurations should depend on initial compaction dynamics: increasing the rate of compression could produce denser equilibrated packings, just like a larger height of free fall, whence a larger initial kinetic energy, increases the density of configurations obtained by rain deposition under gravity Emam et al. (2005). We request the reduced compression rate , defined in (23), not to exceed a prescribed maximum value . The choice of and yields solid fractions and backbone coordination numbers , with a rattler fraction . These data correspond to  kPa (or ). Very similar results are obtained on using a different, but low enough pressure, such as or even  kPa, as remarked in Somfai et al. (2005) (where 2D samples were assembled by oedometric compression), and as indicated in table 3. However, a quasistatic compression from  kPa to or  kPa produces slightly different states at the same pressure. The influence of should disappear in the limit of slow compression, . A practical definition of a (-dependent) limit of loose packing obtained by direct compression, can therefore be proposed as the limit of our D states.

Origin P (kPa) (%) (%)
Table 3: Isotropic states of type D, from direct compression of the granular gas at the indicated pressure (rows marked “gas”), or from gradual, quasistatic compression (rows marked “QS”) of solid samples made at the lowest pressure 1 kPa (or ). Tests of the influence of viscous dissipation parameter and maximum value of reduced strain rate in compression are also made for configurations compressed from a granular gas to 10 kPa.

As reported in table 3, a value of the damping parameter ten times as small as the standard one results in quite similar configuration properties, on compressing a loose granular gas under  kPa, with . So did in fact faster compressions, with , keeping . The data of Table 3 thus suggest that we very nearly achieved the independence on dynamical parameters that is expected in the limit with our choice of control parameters. We note, however, that other possible definitions of a random loose packing, such as the one by Onoda and Liniger Onoda and Liniger (1990) result in different (smaller) solid fractions. Looser arrangements of equal-sized spherical particles can also be stabilized with adhesive contact forces, e.g. on introducing the capillary attractions produced by the menisci formed by a wetting fluid in the interstices between neighboring grains Kohonen et al. (2004).

In addition to packing fraction , coordination number , fraction of rattlers , Table 3 lists the reduced second moment of the normal force distribution, as defined in  (30), the proportion of two-coordinated beads (to be discussed in Section IV.4), , and the average values of ratios (friction mobilization) among contacts carrying normal forces larger and smaller than the average, respectively denoted as and . As a result of some amount of quasistatic compression of the initial assembly, the width of the force distribution decreases, as witnessed by smaller values of in table 3, and so does the mobilization of friction, as measured by and . The effects of compression on the structure and the forces are further studied in paper II Agnolin and Roux (2007a). On comparing numerically simulated loose packings to experiments, it should be recalled that samples are assembled under low pressures in the laboratory: the hydrostatic pressure under a 1 cm thick layer of glass beads is about  kPa. Numerical configurations under higher confining pressures corresponding to mechanical tests in the laboratory (e.g., sound propagation) are more appropriate models if the testing pressure is significantly larger than the initial, assembling pressure – as for the “QS” samples of table 3.

The effects of such proportions of rattlers in granular packings as reported in table 3 have to our knowledge never been studied in detail. It should be emphasized that this relatively large population of rattlers does not jeopardize the global stability of equilibrium configurations, as the stiffness matrix of the force-carrying network is found devoid of floppy modes (apart from harmless, localized ones associated with two-coordinated particles, to be discussed in Section IV.4).

Our D samples should be compared with the simulations reported by Zhang and Makse Zhang and Makse (2005), in which loose sphere packings were also prepared by isotropic compression. Those authors observed, in some cases, lower packing fractions than D values, . Their assembling method is however different: they use a strain-controlled procedure, with a constant compression rate, and then relax the final state at constant volume. In this approach, the pressure reaches very high levels, several orders of magnitude as large as the final value, before samples finally stabilize (Zhang and Makse, 2005, Fig. 3). Zhang and Makse report some dependence of the final state on the compressing rate. Once translated into the dimensionless parameter we have been using here, strain rates used in Zhang and Makse (2005), defined with the typical pressure value  kPa, range between and . The slowest compression reported in Zhang and Makse (2005) is therefore 100 times as fast as the upper limit for we have been enforcing in this work. Viscous forces also differ between the present simulations and those of Ref. Zhang and Makse (2005), in which “global damping” terms are used (i.e., forces opposing the individual motion of particles, rather than relative motions).

Use of low friction coefficients: imperfect lubrication

One way to limit the effects of friction consist in lubricating the grains, as in the experimental study reported in Jia and Mills (2001). If all intergranular friction could be suppressed in the assembling procedure, i.e., for perfect lubrication, then the structure of isotropic packings would be the one denoted as A, studied in Section III. The effect of a small friction coefficient in the contacts while the grain assembly is being compressed can be regarded as a crude, simplified model for imperfect lubrication. We made samples, denoted as B, by compressing the granular gas, just like in the A and D cases, with . In order to approach the limit of slow compression rates better, we started from D configurations, decreased the friction coefficient to , and then requested that while the samples got further compressed to equilibrium under 1 kPa (). (In view of the results in the D case of Section IV.2.1, we do not expect the final B state to be sensitive to damping parameter .) We observed that this small friction coefficient had a notable effect on the final solid fraction, as the value is significantly below the frictionless (A) result, while the coordination number on the active structure is slightly reduced, down to , and the fraction of rattlers raised slightly, to .

Dense, frictional packings obtained by shaking

Another practical strategy to obtain dense configurations is to shake, vibrate or apply repeated “taps” on granular samples Nowak et al. (1998); Philippe and Bideau (2002); Richard et al. (2003). Such procedures involve the introduction of kinetic energy into already quite dense assemblies. In order to investigate their possible effects at a limited computational cost, we avoid the direct simulation of repeated shakes and adopted the following procedure. Starting from the dense A configurations (made without friction and described in Section III), we first apply a homogeneous expansion, multiplying all coordinates by a common factor slightly larger than . With equilibrated A states under confinement level , the chosen value is more than enough to separate all pairs in contact. Then, in order to mimic, in an idealized way, the motion set up by a shaking excitation, the beads are given random velocities (chosen according to a Maxwell distribution), and interact in collisions which preserve kinetic energy, while the volume of the cell is kept constant. This “mixing” stage is simulated with the “hard sphere molecular dynamics” (event-driven) scheme (just like our initial granular fluids are prepared at , as described in Section III.2.2). It is pursued until each particle has had collisions on average. The final preparation stage is a fast compression: velocities are set to zero, particles regain their elastic and dissipative properties (as defined in Section II.1, with friction coefficient , and viscous dissipation, ), the external pressure  kPa is applied via the deformable periodic cell, until a final equilibrium is reached.

The final state is hereafter referred to as C. Quite unsurprisingly, its solid fraction, , stays very close to the RCP value obtained in the A state. However, the coordination number is considerably lower, , which is as small as the value obtained in the loose () D state, while the proportion of rattlers raises to . Remarkably, on comparing states B and C, the latter has a higher density, but a much lower coordination number, and a much higher fraction of rattlers.

We did not thoroughly investigate the influence of parameters and