A Crystallization under shear

Frictionless bead packs have macroscopic friction, but no dilatancy.


The statement of the title is shown by numerical simulation of homogeneously sheared assemblies of frictionless, nearly rigid beads in the quasistatic limit. Results coincide for steady flows at constant shear rate in the limit of small and static approaches, in which packings are equilibrated under growing deviator stresses. The internal friction angle , equal to degrees in simple shear, is independent of average pressure in the rigid limit and stems from the ability of stable frictionless contact networks to form stress-induced anisotropic fabrics. No enduring strain localization is observed. Dissipation at the macroscopic level results from repeated network rearrangements, like the effective friction of a frictionless slider on a bumpy surface. Solid fraction remains equal to the random close packing value in slowly or statically sheared systems. Fluctuations of stresses and volume are observed to regress in the large system limit. Defining the inertial number as , with the grain mass and its diameter, both internal friction coefficient and volume increase as powers of in the quasistatic limit of vanishing , in which all mechanical properties are determined by contact network geometry. The microstructure of the sheared material is characterized with a suitable parametrization of the fabric tensor and measurements of coordination numbers.

45.70.-n, 83.80.Hj, 81.40.Lm, 83.10.Rs

I Introduction

Packings of particles appear in a variety of fields of condensed matter physics and material science, such as granular materials Herrmann et al. (1998); Hinrichsen and Wolf (2004); García Rojo et al. (2005), powders Castellanos (2005), or concentrated, non-colloidal suspensions Stickel and Powell (2005); Ovarlez et al. (2006). Such systems are macroscopically disordered, and share many common features in their rheological behavior. One is a certain shear stress threshold, above which they roughly qualify as a fluid, and below which they might be regarded as solid. In assemblies of particles with purely repulsive force laws, interactions often do not introduce any stress scale, and the threshold only involves some ratio of stress components, whence a behavior often expressed as a friction law. Another basic property shared by many particulate systems is the existence of a specific value of the particle density, above which the material cannot flow. The viscosity of a dense suspension diverges as the solid fraction approaches some value , often regarded Brady (1993) as identical to the random close packing one, ( for identical spherical balls Cumberland and Crawford (1987)). Shearing and volumetric strains are coupled in granular media, which, once densely packed, cannot deform without expanding: this is the dilatancy property, first introduced by Reynolds Reynolds (1885). Once the shear strain reaches a large enough value, granular packs can continuously deform, like perfectly plastic materials, under constant stresses while keeping a constant solid fraction : this state of steady plastic flow does not depend on initial conditions and is known in soil mechanics as the critical state Wood (1990). Friction and dilatancy are coupled in granular materials by the stress-dilatancy relations, as proposed, e.g. by Rowe Wood (1990); Rowe (1962).

It is tempting to try and identify simple, model systems apt to explore the microscopic origin of those broadly defined rheological features. To this end, discrete particle numerical simulation, for granular materials Cundall and Strack (1979); Herrmann et al. (1998); Roux and Chevoir (2005), or suspensions Brady and Bossis (1988), has now become a widespread research tool. Thus friction laws in model granular assemblies in steady shear flows, with some inertial effects, were simulated da Cruz et al. (2005); Hatano (2007), and stress-dilatancy relations were tested Taboada et al. (2006). Many results were obtained on sphere packings Thornton (2000); Suiker and Fleck (2004), which, long investigated in order to characterize their geometry Cumberland and Crawford (1987); Bideau and Hansen (1993), are now studied with complete mechanical models. Thus it has been checked O’Hern et al. (2003); Donev et al. (2005); Agnolin and Roux (2007a) that the random close packing state of monosized spheres is apparently uniquely defined if enduring agitation inducing traces of crystalline order is avoided in the assembling stage. The macroscopic (or internal) friction coefficients, and their relation to micromechanical parameters, including intergranular friction, have been evaluated from numerical simulations Thornton (2000); Suiker and Fleck (2004).

Despite recent advances, some open gaps and unanswered questions can be pointed out in the literature. The accurate and detailed characterization of frictionless systems under isotropic loads O’Hern et al. (2003); Donev et al. (2005), in which static equilibrium states are studied, and few parameters are introduced, contrast with the more general investigations of the behavior of granular systems with intergranular friction Thornton (2000); Suiker and Fleck (2004); Taboada et al. (2006), which are most often carried out by dynamical methods involving inertia effects, and involve quite a few additional parameters. In those latter studies, the limit of frictionless grains is not really treated with the desirable accuracy. Yet, frictionless packings, albeit reported to exhibit singular properties Combe and Roux (2000); Roux and Combe (2002); Tkachenko and Witten (2000), incorporate basic geometric effects that are common to suspensions and dry granular systems, even though they are supplemented by viscous flow effects in the first case, by intergranular friction and inertia in the second case.

In order to clarify issues that have not been settled, the present paper is devoted to a numerical study of frictionless bead packings, subjected to homogeneous shear, and addresses the following questions. Can frictionless packs sustain shear stresses in static equilibrium states as well as in slow, steady flow, and do static and dynamic friction coefficients coincide ? Do fluctuations on measured stresses or strain rates regress in the large system limit ? What can be said about characteristic densities , ,  ? How do classical approaches of dilatancy Reynolds (1885); Goddard and Didwania (1998), and the way it couples to friction Taboada et al. (2006), apply in such a simple case ?

The paper is organized in four main parts. Section II describes the model material and the numerical simulation setup, specifying the boundary condition and initial states used in static and dynamic approaches. Section III reports on the main results about the macroscopic behavior – macroscopic friction and dilatancy – and their dependence on the dimensionless control parameters identified in Section II. Section IV investigates the packing microstructure and the force networks, in connection with macroscopic mechanical properties, with, in particular, a detailed characterization of anisotropy. Section V is a discussion.

Ii Model material and numerical experiments

ii.1 System and interactions

We consider packings of equal-sized spherical beads of diameter and mass , enclosed in a cuboidal simulation cell.

Beads interact in their contacts where only normal forces are transmitted, which are modeled as a sum of an elastic term and a viscous one, as in many numerical studies of granular systems (see e.g., Refs. Thornton (2000); Silbert et al. (2002a); Agnolin and Roux (2007a); Xu and O’Hern (2006)). The elastic force is related to the normal deflection of the contact by the Hertz law Johnson (1985),


where is a notation for , is the Young modulus of the solid material the spherical grains are made of, and its Poisson ratio. Eq. (1) attributes to contacts a variable spring constant .

The viscous normal force opposes the relative normal velocity of contacting beads, and is chosen as


with a constant coefficient . The same form of the viscous force was used in Somfai et al. (2005); Agnolin and Roux (2007a). Although (2) is devoid of a physical justification, some kind of dissipation is required (a granular material is not a conservative system), and consequently, the influence of on the simulation results has to be carefully assessed. One attractive feature of the force law (1) and (2) is the resulting velocity-independent coefficient of restitution in binary collisions. Most simulations reported here were done with values such that is close to zero.

Particle rotations play no role and are ignored, as frictionless spherical objects behave like point particles interacting with central forces.

The equations of motion for the particles, given by Newton’s law, as in all molecular dynamics (MD) methods, are to be numerically solved with standard time discretization schemes Allen and Tildesley (1987). The time step used in the computations is a small fraction of the characteristic period of oscillations for the stiffest contact.

ii.2 Boundary conditions, stress and strain control

We use different simulation procedures in which some strain, or strain rate, and stress components are externally imposed to the system.

In order to avoid wall effects and to determine easily the intrinsic constitutive laws that apply in the large system limit, the simulation cell has periodic boundary conditions. The edges of the cell have lengths along the three orthogonal axes of coordinates. Unlike the cell, the material may undergo some shear strain, imposed with the Lees-Edwards procedure Allen and Tildesley (1987). Adding this possibility to the potential shrinking deformations along the three axes of coordinates, four independent strain components are considered in the different simulation steps and methods we are using in this work. The procedures defined below consist in choosing to fix some of them to zero or to a constant value while prescribing the values of stresses conjugate to the others. Table 1 recapitulates those choices for the three different simulation procedures.

Initial assembling process: procedure O

In a preliminary step, the system is first prepared by isotropic compression of a loose “gas” of particles. The corresponding procedure, denoted as “O” (like “origin”), is the one applied in Agnolin and Roux (2007a) to prepare isotropic packings. Global shear strain is kept equal to zero, while the system shrinks along all three directions, until a mechanical equilibrium state is reached for which all three diagonal components of the Cauchy stress tensor Bathurst and Rothenburg (1990); Christoffersen et al. (1981) are equal to a set pressure value . The system is deemed equilibrated when all forces compensate to zero, with a tolerance set to on each particle, and each diagonal stress component is equal to with a relative error smaller than , while the kinetic energy per particle does not exceed . Those isotropic equilibrated configurations are the “random close packing states”, as studied in O’Hern et al. (2003); Donev et al. (2005); Agnolin and Roux (2007a).

Controlled shear rate: procedure D

Initial configurations produced with method O may then be subjected to a simple shear deformation, in which a macroscopic motion along direction 1 is set up, with velocity gradient, on average, along direction 2 (by the Lees-Edwards procedure), while and are fixed. is allowed to fluctuate in order to maintain equal to on average (with very small fluctuations). The macroscopic shear rate is denoted as . This defines procedure “D” (for dynamically sheared). It was implemented in a very similar way in Rognon et al. (2006). One then records the time-averaged shear stress , as well as the sample volume. It is important to note that Lees-Edwards boundary conditions are fully compatible with either a linear velocity profile or very heterogeneous strain fields. With this procedure, shear strain is set equal to the ratio of the offset along axis of the neighbor copy of the simulation cell in the direction, to the length . Consequently, due to fluctuations in , the time derivative of is not strictly equal to at all times.

Static approach, controlled shear stress: procedure S

In the limit of small , results of procedure D simulations should be comparable to static computations, in which the system equilibrates under an externally imposed shear stress. To compare static and dynamic measurements (possible differences between “static” and “dynamic” friction coefficients or threshold shear stresses in similar systems are mentioned in Hatano (2007), and discussed in Da Cruz et al. (2002); Xu and O’Hern (2006)), we also implemented a completely stress-controlled, quasistatic procedure, denoted as “S” (for static). In procedure S, increasing values of shear stress are stepwise applied, by increments , to the initially isotropic configurations obtained with procedure O, while the prescribed value of all three diagonal components is the initial pressure . , unlike in procedure D, is not constant. It satisfies a dynamical equation designed to impose a prescribed value to . For each value of , one waits until a satisfactory equilibrium state is reached (with the same tolerance levels as in procedure O). The calculation is stopped when the packing does not equilibrate with the current value of after MD time steps. The largest value for which an equilibrium state was obtained is kept as an estimate of the shear stress threshold for onset of flow.

Stress/strain control Procedure O Procedure D Procedure S
Table 1: Choice of imposed stresses or strain rates in the three simulation procedures O, D, and S.

ii.3 Dimensional analysis, state parameters and geometric limit

Assuming homogeneous steady states are observed in large enough samples, with shear rate and normal stress , then, by dimensional analysis Roux and Chevoir (2005); da Cruz et al. (2005); Agnolin and Roux (2007a) all dimensionless state variables, such as solid fraction or average stress ratio only depend on three dimensionless parameters.

The first one, the inertial number, , characterizes the importance of inertial effects in dense granular flows GDR MiDi (2004); da Cruz et al. (2005); Hatano (2007) and plays a central role for these systems Jop et al. (2006); Agnolin and Roux (2007b). The quasistatic limit is the limit of , which we will systematically explore.

The importance of contact deformation is characterized by the second dimensionless parameter, a stiffness number which we define, like in Agnolin and Roux (2007a), as . is such that the typical contact deflection under pressure is proportional with a prefactor close to 1 Agnolin and Roux (2007a). In order to enable comparison of macroscopic elastic properties with experimental results, we set  GPa and (glass elastic constants). The pressure levels chosen,  kPa and  kPa, then respectively correspond to and . These two values of were reported to be large enough for the limit of rigid grains, i.e., of , to be approached with good accuracy in the case of static packings Agnolin and Roux (2007a).

Finally, the third dimensionless parameter is the level of viscous damping , which appears in a viscous force and should not play a major role in the quasistatic limit.

Table 2 sums up the values of dimensionless control parameters used in the present numerical study.

Table 2: Range of dimensionless parameters used in this study.

We should investigate the relations between global intensive variables, such as stresses, density, strain rate, in the limit of large samples, i.e., of . It is expected that for large enough samples the material state in shear flow will not depend on the specificities of boundary conditions, or on whether shear stress or strain rate is controlled. This requires the investigation of possible size effects and the study of the regression of fluctuations for global variables. Measured state variables should also be uniform in space – and thus one needs to check for possible shear localization. If dimensionless variables such as stress ratios or density are well behaved in the triple limit of (thermodynamic limit), (quasistatic limit) and (rigid limit), then the observed inner states and mechanical behavior of the packings only depend on their geometric properties – hence the name macroscopic geometric limit we adopted for such a situation. One of the major goals of the present study is the investigation of material properties in this limit.

Finally, as a practical application of the dimensional analysis of simulation parameters, let us note that the computational cost, expressed as a number of MD integration steps needed to reach a given shear strain , is proportional to .

Iii Global variables and macroscopic behavior

Our global observations and measurements are reported in this section. Conditions for proper observations of the intrinsic behavior of the material subjected to procedure D (shear-rate-controlled numerical experiments) are checked for in Section III.1, in which various qualitative aspects of the material state in shear flow are discussed. Attention is then focused on macroscopic friction (Sec. III.2) and dilatancy (Sec. III.3) properties of the material, which are more thoroughly and quantitatively investigated. Finally, the results obtained with procedure D at low inertial numbers are compared to those of the static approach, procedure S, in Section III.4. Section III.5 discusses the essential results and their connections with the literature on granular materials.

Figure 1: (Color online) (left axis, in black) and (right axis, in red) as functions of strain . Note that the left and right scales are different. Time series obtained with , , and .

iii.1 Material state in slow shear flow: qualitative aspects

With procedure D, we investigate steady states, and time series are collected for averaging. We are interested in intrinsic constitutive laws, as measured on averaging over the whole sample. It is therefore necessary to check for both invariance in time and homogeneity, in the statistical sense. We should also assess the control of constant stress , and discuss the values of other stress components.

Steady state flows and stress measurements

Fig. 1 displays the evolution of two components of the stress tensor, and , with strain . It shows that is well controlled since it was requested to stay equal to in this numerical experiment. The evolution of stress , from the initial, isotropically confined state, witnesses the existence of an initial transient, which has virtually ended at in that case. The steady state part of the time series starts for values of strain that depend on the inertial parameter, of order for the smallest values, (), increasing typically to about for and to several units for . Unlike in dense systems with intergranular friction Thornton (2000); Suiker and Fleck (2004); Taboada et al. (2006), for which deviator stresses, starting from isotropically compressed initial states, go through a peak before approaching a plateau value at large strain, the shear stress in frictionless bead packs appears to grow monotonically, as a function of strain, toward its steady state value. Another notable feature of the shear stress as a function of time is the importance of fluctuations, which often exceed 30% of the mean value on the example of Fig. 1, in a sample of 4000 beads. A proper evaluation of average shear stresses thus requires careful statistical approaches and error estimates.

As a practical criterion to detect the end of the initial transient regime, we request that a small set of basic measured quantities do not exhibit any visible trend. Specifically, shear stress , volume fraction and coordination number should all fluctuate about their mean value in a stationary manner, as well as the kinetic energy per particle, , associated with velocity fluctuations. The latter is defined as


measures the instantaneous discrepancy between the actual flow generated by the Lees-Edwards boundary condition in the granular material and the affine velocity field in a homogeneous continuum in shear flow.

Unlike , lengths and are constrained to remain constant in procedure D, so that and may vary during the simulation. For , we observed that time averages of and differed from the initial hydrostatic pressure by less than . This difference becomes even smaller for smaller inertial numbers: for , relative differences and respectively reduce to and . Those values decrease down to and for , and to and for . Although apparently not equal to zero, even in the quasistatic limit, those stress components are very small, and, consequently, will not be studied in the sequel. Sec. III.2, instead, focuses on accurate determinations of shear stress .

For a given number of particles, the relative fluctuations of the instantaneous value of , and (i.e., the ratio of their quadratic average to the mean value) seem to be independent of . The average values of , on the other hand, as compared to the kinetic energy of the macroscopic field, which is proportional to , increases as decreases. Fig. 2 is a plot of versus , showing that this ratio approximately diverges as in the limit of . This agrees with measurements made in 2D simulations of shear flows: the same behavior is reported in Ref. da Cruz et al. (2005), and an interpretation was suggested, to which we shall return in Section III.5.

Figure 2: Kinetic energy associated with velocity fluctuations, as defined in (3), normalized by , versus , in simulations with 4000 beads, for and .

These observations suggest that in the quasistatic limit one has increasingly inhomogeneous instantaneous velocity gradient fields, which we now investigate.

Instantaneous velocity profiles

Instantaneous velocity profiles recorded at different random times for different values of are plotted in Fig. 3.

Figure 3: (Color online) Two velocity profiles at randomly chosen times, for (top), (middle), (bottom). , and .

Profiles are obtained on averaging particle velocities over slices cut alongside in the simulation cell (particles overlapping slice boundaries contribute to several different averages). Inertial number has an important effect on the granular flow. As shown in Fig. 3, instantaneous velocity profiles for are linear, whereas shear bands may appear for , as in the profile marked “L” (for “localized”) on the bottom plot of Fig. 3. The transition between these two regimes seems to be gradual, with profiles in the middle part of Fig. 3, corresponding to , exhibiting somewhat intermediate features.

Localization occurs here in the bulk of the material since the system is not enclosed between walls. Localization is thus an intrisic property of the studied material, which spontaneously appears for small values of  Aharonov and Sparks (2002).

At first glance, it seems that the erratic behavior of the velocity profiles in the quasistatic limit may seriously jeopardize the interest of the results obtained on averaging over the whole simulation size and would demand specific analysis, distinguishing between material states within and outside shear bands. However, localization patterns are not persistent, and linear velocity profiles are recovered by time averaging, even in the limit, which means that on average, the flow is homogeneous. Figure 4 shows the gradual fading out of strain rate localization, after a strain interval of order . Shear bands thus randomly appear, move and disappear. Such a behavior is witnessed by larger relative fluctuations of as decreases. We did not carry out a detailed study of the lifetime and dynamics of nonpersistent shear bands, as the statistical homogeneity of the system in steady state shear justifies an analysis of global behavior based on averages over space and time.

Figure 4: (Color online) Velocity profile after shear strain intervals equal to 0.004 (red dotted curve), 0.02 (blue dashed curve) and 0.1 (black solid curve), following the instant corresponding to the localized profile marked L in Fig. 3, bottom graph.

iii.2 Macroscopic friction coefficient

For D simulations, the macroscopic friction coefficient, which we denote as , is set equal to the time average – in the steady state – of the ratio of the shear stress to the normal stress (we use a convention where compressive stress components are positive)


The simulations produce raw data in the form of time series. The steady part of the time series is isolated as explained in Sec. III.1 and can then be easily computed. To estimate the statistical uncertainty on the measurement of averages over finite time series, we use the “blocking” (or “renormalization group”) technique presented in Flyvbjerg and Petersen (1989). This yields error bars on measurements of averages in finite systems which should not be confused with the quadratic average of fluctuations of the observable quantity. In practice, due to intrinsic long-lasting correlations in the system, we observed that quite long runs were necessary. In some cases with , up to simulation time steps (corresponding to a deformation ) were necessary for a correct evaluation of the uncertainty on .

In the present case, we could check that the time series of all observable quantities recorded in distinct samples differing only by their initial state were statistically identical in steady state with high accuracy, as expected from critical state theory Wood (1990); Radjaï and Roux (2004); Rothenburg and Kruyt (2004); Kruyt and Rothenburg (2006).

, as estimated from time series in type D simulations, may depend on the three dimensionless numbers introduced in Sec. II.3, as well as on the number of particles. This dependence is investigated in the following paragraphs.

Effect of

Among the three dimensionless parameters governing the behavior of the system, the inertial number has the strongest effect on . Fig. 5 plots as a function of . It shows that is an increasing function of . This dependence of the macroscopic friction coefficient on the inertial number is similar to the ones reported in the literature, as obtained by both simulations and experiments GDR MiDi (2004); da Cruz et al. (2005); Hatano (2007), although most published results pertain to granular systems with friction in the contacts. Here approaches a finite nonzero value in the quasistatic limit of , despite the absence of friction at intergranular contacts. coincides with the internal friction coefficient of the material in its critical state.

Note the accuracy of the displayed curve: statistical uncertainties measured with the blocking method are comprised between and and are thus invisible on the graph.

Figure 5: Macroscopic friction vs. inertial number for stiffness parameter , damping parameter and number of beads . The solid line is Eq. (6) with the parameters of Table 4 (no visible difference on using best parameters with or ).

Effect of

Near the rigid limit, the macroscopic behavior should reflect the absence of stress scale in the contact law: stress ratios and derived quantities such as the macroscopic friction coefficient are expected to be independent of the average stress. The friction coefficient hardly changes between the two values of used in our simulations, indicating that the rigid limit is accurately approached. Those simulations were carried out with for , and the relative variation on is less than throughout this range of inertial parameter on varying the stiffness parameter from to . Thus we decided to gather the values obtained for the macroscopic friction coefficient with and in Fig. 5, because the uncertainty on the macroscopic, geometric limit of to be estimated will eventually exceed this small difference.

Effect of

The viscous damping term is indispensable in the model, as the only source of dissipation, but its magnitude should be irrelevant in the quasistatic limit. Consequently, the influence of on our results had to be assessed and we performed simulations for different values of with (this value corresponds to a restitution coefficient ), (), () and (). Our results show that for , the maximal relative variation of on the macroscopic friction is less than . Far from the quasistatic regime, the influence of is no more negligible: the relative variation of is greater than on changing when .

Effect of

The influence of the sample size on the average of the apparent friction coefficient , was investigated on comparing results for three different numbers of particles: , and . Results are listed in Table 3. We also recorded the standard deviations, denoted as , and the average of the top percentile of the instantaneous values, denoted as . Let us recall that we are dealing here with the fluctuations of the time series, which differ from the statistical uncertainties on the average values.

500 0.1169 0.3100 0.2188 0.6367 0.0022 0.6403
1372 0.1101 0.1907 0.1609 0.6380 0.0015 0.6408
4000 0.1090 0.1245 0.1431 0.6387 0.0008 0.6404
5002 0.1432 1.263 0.8378 0.6738 0.0178 0.7148
1372 0.1209 0.1689 0.1779 0.6365 0.0016 0.6390
4000 0.1197 0.1002 0.1519 0.6368 0.0010 0.6388
500 0.1473 0.2091 0.2275 0.6316 0.0027 0.6360
1372 0.1457 0.1293 0.1966 0.6322 0.0016 0.6346
4000 0.1458 0.0764 0.1752 0.6323 0.0009 0.6338
500 0.2112 0.2045 0.3317 0.6193 0.0025 0.6234
1372 0.2123 0.1197 0.2802 0.6197 0.0015 0.6223
4000 0.2125 0.0694 0.2517 0.6200 0.0009 0.6215
Table 3: Influence of sample size on macroscopic friction and volume fraction for different values of inertial number , with and . Superscripts “” denote the average of the top percentile values in the steady state part of the time series.

The effect of the sample size on the macroscopic friction is unnoticeable for and since the difference between the friction coefficients pertaining to these two sizes is less than the statistical uncertainty marring the accuracy on . However, the impact of on cannot be neglected in the quasistatic limit for a system of beads. These results show that some minimum number of beads, of order about 1000, is required to approach the thermodynamic limit with an acceptable accuracy.

The data of Table 3 also witness the regression of fluctuations of stress ratio in the steady state in the large system limit. The results are compatible with the classical form for the decrease of fluctuations of collective variables, viz. . Specifically, for , and , a fit of the data to the following form:


has good statistical admissibility. This result is shown in Fig. 6 in graphical form (two additional sizes and were also simulated).

Figure 6: as a function of for , and . The solid line equation is relation (5).

Therefore, we expect the steady state stress-strain curves such as the ones of Fig. 1, however noisy for the sample sizes simulated, to become smooth in the large system limit.

Approach to the macroscopic geometric limit

According to the previous parametric study, the geometric limit can be confidently studied as the limit of on samples of beads with and .

should be a function of the sole inertial number in very good approximation for sufficiently small values of . In the absence of any scale, we tried to fit the data by a power law function (see Fig. 5) of the form


As stated above, this fit is not expected to be accurate for large values and we therefore restricted ourselves to fit the data points with . Parameters (the geometric macroscopic friction coefficient), and were separately estimated for and (keeping and ) and results are shown in Table 4.

Table 4: Best fit parameters for Eq. (6) and the data obtained with , for and .

The value of the geometric macroscopic friction angle corresponding to is (for )


Quite similar values are also reported with two-dimensional packings of polydisperse disks by Taboada et al. Taboada et al. (2006), whose estimate of the macroscopic friction angle lies between and for frictionless grains, and by Da Cruz et al. da Cruz et al. (2005), who obtained in shear flow simulations for small parameters. Hatano Hatano (2007) recently performed 3D numerical simulations on polydisperse assemblies of about 10000 spherical beads, for different intergranular friction coefficients . The reported value of the macroscopic friction coefficient in the quasistatic limit is for , apparently lower than our result. It should be recalled however that Hatano’s work was motivated by applications to granular materials under high confining stresses within geological fault zones, and that consequently simulations were carried out with lower stiffness levels (, , and ) than in the present study. Moreover, the lower range of parameters, below , was only investigated with the lower values. Hatano used the same form as Eq. (6) to fit his data, but his estimate differs from ours (see Table 4). Although some effect of the polydispersity is possible, we also attribute this discrepancy to some non-negligible influence of in Hatano’s simulations Hatano (2007). Only the simulations with in Hatano (2007) can be expected to approach the rigid limit accurately. For this stiffness level, Hatano’s data points are available for and are in very good agreement with ours (e.g., for ).

iii.3 Dilatancy and steady-state density

Dilatancy under shear is a basic property of granular materials in quasistatic deformation Reynolds (1885); Wood (1990); Rowe (1962); Taboada et al. (2006), when dense samples are subjected to increasing deviator stresses. The steady-state density is mainly sensitive to if is large enough da Cruz et al. (2005). The small behavior of frictionless bead assemblies is investigated here with greater accuracy than in previous studies.

We could check that, just like the macroscopic friction coefficient, the steady state time average of the volume fraction, , is an intrinsic property of the material, independent of its initial preparation. Next, we investigate its dependence on the three dimensionless parameters , and and on the number of particles .

Effect of , and

Once again, numerical experiments demonstrate that among the three dimensionless numbers governing the behavior of the system, the inertial number has the most important effect on . Fig. 7 shows the influence of on . We observe that decreases for increasing , as previously reported da Cruz et al. (2005); Hatano (2007). It starts from a value in the quasistatic limit and the system expands as increases. Statistical uncertainties on measured thanks to the blocking method are comprised between and and are thus invisible on the figure.

Figure 7: (Color online) Volume fraction as a function of inertial number (for , ), for both stiffness levels (blue squares) and (red triangles) The solid lines are given by Eq. (11) with the parameters of Table 5.

is very close to  O’Hern et al. (2003); Agnolin and Roux (2007a), which coincides (up to small corrections due to the finite contact stiffness) with the initial volume fraction , right after the samples are assembled with procedure O. For and we have , and for (averages and standard deviations are evaluated on five samples). The system studied thus appears to be devoid of dilatancy in the quasistatic limit. Whether should be regarded as equal to at the macroscopic level will be discussed after the possible influence of on the average densities is investigated.

Stiffness parameter typically induces a relative increase of the volume fraction of roughly when it changes from to , whatever the value of – a small effect, yet distinguishable from statistical uncertainties. Such a density increase is of course expected, as larger contact deflections due to larger stresses or a softer material decrease the sample volume. Simulations with (), (), () and () for a wide range of inertial numbers have also been run. The influence of on is important for large : for , the relative variation of with can reach . However, this effect, as expected, gradually vanishes as the quasistatic limit is approached, and for the relative variation of with is less than .

Effect of

According to Table 3, very slightly varies with the number of particles, like in static, isotropic systems O’Hern et al. (2003); Agnolin and Roux (2007a). The following fit, based on the measurements for the smallest available value of , i.e., for , may be used:


with the following parameters:


As with the macroscopic friction coefficient , we could check for the regression of fluctuations of the volume fraction for increasing . For the same set simulations with , and as in Sec. III.2.4, we observe that the decrease of density fluctuations with increasing can be fitted by the following relation:


as shown graphically in Fig. 8,

Figure 8: as a function of for the same time series as in Fig. 6, fitted with Eq. (10) (solid line).

whence a well-defined large system limit for .

Approach to the macroscopic geometric limit

Volume fraction can therefore be modeled as a function of near the quasistatic limit, say for , with small corrections to account for the influence of and . It can be regarded as independent of (at least for ). The fit form used is


For and , the best fit values of the parameters of Eq. (11) are given in Table 5.

Table 5: Best fit parameters for Eq. (11) and the data obtained with , for and .

To evaluate the macroscopic value in the double limit of and , it is reasonable to assume that the small corrections to that result from the finite value of and from the nonvanishing value of are additive. The use of Eqs. (8)-(9) to evaluate the finite correction to the value of the quasistatic density, as obtained on fitting Eq. (11) to the results with , yields, for :


The increases of , from its value in the rigid limit, due to the finite stiffness is of order (see (Agnolin and Roux, 2007a, Eq. 31)) and is smaller than the statistical uncertainty in (12). The value of given in (12) is thus our best estimate, from D-type simulations, of the solid fraction of sheared sphere packings in the macroscopic geometric limit.

iii.4 Static behavior

We now compare the results of Sections III.2 and III.3 for steady shear-rate-controlled simulations (procedure D) with those obtained through static shear numerical experiments (procedure S).

Friction coefficient

The static macroscopic friction coefficient is defined in procedure S as


where denotes the maximum shear stress which the system has been able to sustain in mechanical equilibrium, and the confining pressure. Static microscopic friction coefficients for different sample sizes are displayed in Table 6.

256 4 0.246 0.022
500 4 0.210 0.007
1372 6 0.169 0.004
2048 6 0.154 0.004
2916 6 0.145 0.007
4000 10 0.136 0.007
8788 6 0.122 0.005
Table 6: Average and standard deviation of the static friction coefficient obtained in S-type simulations, over samples of grains for different . Data corresponding to both values of (with samples each) are aggregated.

Values of are larger than the dynamical value obtained in D simulations in the quasistatic limit. As shown by Tab. 6, is size dependent, unlike (for ). Analogous observations were reported in Xu and O’Hern (2006) for two-dimensional systems of frictionless disks: in the limit of vanishing shear rates, the shear stress reaches its large system limit with only several hundreds of beads, whereas the minimum shear stress required to maintain a long lasting steady shear flow exceeds the previous one and is more sensitive to .

Fig. 9 shows the influence of system size on (discarding the smallest values). The data are correctly fitted by the following relation




The related angle of friction is . This is consistent with the equality, in the thermodynamic limit (), of the dynamical and static macroscopic friction coefficients (see Eq. 7).

Figure 9: (Color online) Size dependence of . denotes the number of particles in the system. The solid line is the fit of Eqs. (14)-(15). Crosses are the top percentile values extracted from the time series of obtained in procedure D, as listed in Table 3. The hashed region represents the estimate, from D simulations, of with its error bar (Table 4). The (blue) dot with an error bar on the left axis is the static estimate, .

The influence of is very small and negligible in comparison to the effect of the system size, and we therefore averaged over systems with both stiffness levels and .

With the smallest system size simulated, , we observed that some of the samples, once submitted to shear stresses, acquired a strongly ordered, crystalline structure, to be discussed in Appendix A.


Static shear simulations with procedure S support the observation made in Sec. III.3 that the frictionless model material studied is devoid of dilatancy in the quasistatic limit. As shown in Fig. 10, which represents as a function of the macroscopic stress ratio imposed to the material in different samples with , the volume fraction hardly evolves with the stress deviator as it is increased towards its maximum value. However, whatever the initial state of the system, it experiences a slight compaction at the beginning of the shear and a small decompaction near the failure limit, but we have no convincing explanation for this phenomena. The evolution of is somewhat erratic (as in previous studies on 2D rigid, frictionless disk assemblies Combe and Roux (2000); Roux and Combe (2002)) and the density change between the isotropic initial state and the one supporting the maximum shear stress is equal to zero, within statistical uncertainties.

Figure 10: Variation of the volume fraction with the static shear imposed to five different samples of 4000 beads with . Each curve stops at a given value of : this is the greatest value for which the packing has managed to reach mechanical equilibrium.

Similarly to the values of measured in D simulations, solid fraction in static packings under maximum shear stress is slightly dependent on sample size, with a negative finite-size correction to the macroscopic value. On fitting a variation proportional to one gets, for ,




values for are slightly larger, by about .

Comparing this estimate of with the result for given in (12), we conclude that static and dynamic solid fractions in quasistatic shear are identical, within statistical uncertainties. Disregarding the very small correction due to the finite value of (equal to about on applying the formula given in (Agnolin and Roux, 2007a, Eq. 31)), this means that, just like for , the values of solid fraction in the macroscopic, geometric limit coincide in strain rate controlled and in shear stress controlled approaches.

As to the value of the solid fraction in the initial isotropic state, a similar evaluation of size effects yields (using the samples of Table 6 with and ):




As announced, this is the random close packing value O’Hern et al. (2003); Agnolin and Roux (2007a). Results (19), (17) and (12) are compatible, and we thus conclude that the system is devoid of dilatancy under shear in the macroscopic geometric limit.

iii.5 Discussion

We briefly review and comment here the essential results on the macroscopic behavior of the material under simple shear, and compare them to other available results on similar systems.

Internal friction and the macroscopic geometric limit

Whether assemblies of frictionless grains have a well-defined, finite internal friction coefficient has sometimes appeared as a debatable issue, although some previously cited works da Cruz et al. (2005); Hatano (2007); Oger et al. (1998) relying on numerical simulations of slow shear flows in steady state agree with our positive conclusion. A proper evaluation of in the macroscopic geometric limit requires more care than corresponding measurements in granular assemblies with friction. This is due to the importance of fluctuations, as apparent on Fig. 1. In D-type simulations, it is also necessary to explore a range of very small inertial numbers to accurately evaluate the quasistatic friction coefficient, as apparent in Fig. 5. As an example, for , quite a small value, the macroscopic friction coefficient already exceeds its quasistatic limit by 25%.

Our estimate for is confirmed by static simulations, once the results are suitably extrapolated to the limit of large systems. One may understand this size effect on S results as follows. The friction coefficient evaluated in D simulations is an average over time series with large fluctuations. However, the system remains close to mechanical equilibrium at any time. Assuming it is possible to find an equilibrated configuration very close to all dynamically explored states, however large the instantaneous value of the shear stress, the S procedure would be able to find statically supported shear stress values as large as the maximum of in D time series. Although clearly oversimplifying the evolution of the system in configuration space, this explanation appears to be correct at least on correlating -dependent maximum static shear stress levels to fluctuations in slow shear flows: the dependence of , as plotted in Fig. 9, is paralleled by that of the typical largest values of (top percentile) in D simulations.

Static and dynamic values of shear stress thresholds for flow are also observed to coincide in the fixed density simulation results of Xu and O’Hern Xu and O’Hern (2006), obtained on 2D packings of frictionless disks, with a similar excess of the static value that vanishes as .

When non negligible inertial effects are present, we observe that the increase of internal friction with inertial number is the dominant feature of the material behavior (the effect of stiffness level is smaller by orders of magnitude), in qualitative agreement with many other results on frictional and frictionless grains da Cruz et al. (2005); Hatano (2007).

Absence of dilatancy

Our results also agree in static shear and in steady state constant shear rate flows for the average volume fraction , which stays equal to its value in the initial, isotropically confined configuration in the macroscopic geometric limit. Our data show that, within statistical uncertainties (i.e., about ) the critical value of is equal to in packings of frictionless spherical beads.

The material studied is thus devoid of dilatancy. Interestingly, this contradicts the simple pictures of the origins of dilatancy which have been proposed since the introduction of this property by Reynolds Reynolds (1885), based on the distortion of simple assemblies of a small number of contacting spheres (like, e.g., a regular tetrahedron) Goddard and Didwania (1998).

The absence of dilatancy in the quasistatic limit is also at odds with the classical ideas on the relation between dilatancy and internal friction, according to which macroscopic friction stems from two microscopic origins, intergranular friction and dilatancy, with an additive combination of relevant angles Rowe (1962); Goddard and Didwania (1998). Ref. Taboada et al. (2006) adds another component to macroscopic friction, due to intergranular collisions as a source of dissipation, and therefore accounts for the internal friction of frictionless grains. Thus is the internal friction angle that we measure in the geometric limit. Ref Taboada et al. (2006), although only incidentally dealing with frictionless materials, nevertheless appears to predict a positive dilatancy in that case, which our results do not confirm. Similarly, a recent study published by Kruyt and Rothenburg Kruyt and Rothenburg (2006), which also deals with 2D disk assemblies, predicts a non-vanishing dilatancy when intergranular friction coefficient approaches zero. Ref. Kruyt and Rothenburg (2006), similarly to Ref. Taboada et al. (2006), discusses stress-dilatancy relations, and finds a linear variation of the dilatancy ratio with the difference between peak and steady-state macroscopic friction. In contradiction with our data, it attributes a positive value to both quantities as the friction coefficient approaches , while its estimate for is significantly larger than our (3D) one, or than the (2D) one of Ref. Taboada et al. (2006). (Note that the maximum deviator to mean stress ratio, as defined in Kruyt and Rothenburg (2006); Radjaï and Roux (2004), is ). The frictionless case was not directly simulated in this work. Some rapid variations of macroscopic friction and dilatancy angles near the singular limit of might be overlooked.

In our simulations, instantaneous fluctuating shear stress and volume fraction, however, appear to be correlated, suggesting some stress-dilatancy coupling at the level of short-lived, transient and rearranging structures, which disappears on taking time averages.

A toy model

Since some of our results on the macroscopic behavior of frictionless bead packs might seem counter-intuitive, we designed a simple model in which similar basic ingredients (geometric constraints defining isolated equilibrium positions, inertia, viscous dissipation) produce an analogous behavior in a suitably defined “macroscopic geometric limit”. In both cases, the microscopic motion is a succession of arrested dynamical phases, alternating with approaches to transient equilibria. The toy model simply provides suggestive analogies, it should not be regarded as a real physical explanation for the macroscopic behavior of the granular system.

We consider a single object of mass , subject to its weight , pushed along a rough horizontal surface. For simplicity the model is two-dimensional, with only one horizontal coordinate, , and the surface profile , along vertical coordinate , is periodic with wavelength , as depicted in Fig. 11. The mobile object is driven either by a constant horizontal force , or by a piston with constant horizontal velocity . Both contacts are rigid and unilateral, so that the mobile object might move faster than the piston if accelerated downhill by gravity, or occasionally take off from the surface.

Figure 11: The model of the slider on a rough surface. (a) Case of a sinusoidal profile. Forces at point S are drawn as vectors. (b) Profile for which .

Force is the analog of shear stress in the granular material, and that of , while horizontal and vertical displacements respectively correspond to shear strain and volume increase. A viscous force opposes the tangential motion along the surface, so that for the slider stabilizes at some local minimum of profile , like point O on the figure. Such minima are analogous to equilibrium states under isotropic pressure.

Let us first discuss the static experiment, in which the mobile object, initially in equilibrium in under , is subjected to a growing horizontal force . It equilibrates where the tangent direction to the substrate is orthogonal to the applied force, . It has first to move upwards, hence some dilatancy. The maximum value of is the static effective friction coefficient of the object on the surface, equal to the maximum slope of profile . It is reached at point S on the figure. Effective static friction angle is the maximum angle between the reaction of the substrate, force on Fig. 11, and the vertical direction. As the quasistatic motion from O to S follows the surface, dilatancy , defined as the ratio of vertical to horizontal coordinates of the velocity (corresponding to ratio in the sheared granular material), is also identical to the maximum profile slope. Dilatancy and friction angles coincide: . If a nonzero friction coefficient is introduced in the contact between the mobile object and the substrate, then reaction at point (see Fig. 11) may form an angle with normal direction , so that the effective static friction angle is – a classical form of the stress-dilatancy relation Wood (1990); Taboada et al. (2006).

In the velocity-controlled case, the dynamic friction coefficient is conveniently evaluated from the dissipation of energy. In the limit of small velocity, the mobile object pulls ahead of the velocity-controlled driving piston at each maximum of . Its subsequent downhill sliding is accelerated by gravity, but it is prevented by viscous dissipation to pass the next maximum, and ends up at the bottom of the valley, where it is later picked up by the slow piston, to be pushed up the next ascending slope. In this scenario the dissipated energy per wavelength is the potential energy loss in a fall over height . Hence an effective friction coefficient . This result is, remarkably, independent of the viscous damping coefficient, just like the macroscopic friction of the granular material in slow shear flow. As the properties of the system only depend then on substrate geometry, the limit of slow imposed velocity is the geometric limit.

The macroscopic limit can be defined as , where is the length scale on which the effective properties of the slider are studied. Consequently, its vertical motion, on the scale of microscopic asperities of the surface, becomes irrelevant, and one observes effective macroscopic friction without dilatancy. Models for dilatancy Goddard and Didwania (1998) apparently focus on microscopic phases of the motion in which the slider rises up the slope, but ignore the equally important ones in which it falls down.

is the average slope of the ascending part of the profile, multiplied by the fraction of length for which is increasing. It is in general smaller than , which is the maximum slope. Thus, for a sinusoidal profile, , as represented on Fig. 11(a), one has , while . In order for both friction coefficients to coincide, function should be as shown on Fig. 11(b), with ascending parts of constant slope, followed by vertical drops.

The particular profile shape of Fig. 11(b) can be argued to be appropriate for the analogy with the bulk material. As long as the contact with the substrate is maintained, the configuration might be an equilibrium position for some (possibly negative) value of . In the analogy with the granular material, the contact network might balance the external load for some value of the applied stress components. The free fall, on the other hand, is the analog of a network rearrangement, during which applied loads cannot be supported because of the missing contacts. In the granular material (as explicitly shown in Combe and Roux (2000)) intervals of stress components for which a given contact network is stable shrink to zero in the large system limit. This corresponds for the toy model to a constant slope of the rising parts of profile (defining a unique possible value of in equilibrium).

Finally, the velocity-controlled sliding of the object on the profile of Fig. 11(b) also provides an interpretation of inertial number  GDR MiDi (2004), and of the behavior of the kinetic energy da Cruz et al. (2005). The motion involves two characteristic times: the duration of the rising phase, in which the object is in contact with the piston and moves with horizontal velocity , ; and the duration of the free fall, . Their ratio, , is the analog of number , as readily checked on replacing distance by some length of order (a typical interstice between neighboring grains to be closed for a new contact to appear), by , and force by , which is the order of magnitude of unbalanced forces on the grains in the dynamical phases of motion. The free fall phases of the motion explain why the kinetic energy is, on average, much larger than the scale associated with the macroscopic motion. More precisely, the time average of the kinetic energy associated with velocity fluctuations is of order (for ), whence (discarding constant factor ) the behavior shown in Fig. 2, .

Iv Microstructure and force networks

Our specific emphasis on the geometric limit of the macroscopic mechanical behavior of frictionless bead packings calls for an analysis of the geometry of sheared configurations, the first motivation of which is to explain the microscopic origins of macroscopic friction. Ultimately, a model should be sought which, unlike the analogous one of Sec. III.5.3, would explicitly and quantitatively describe the mechanisms, involving instabilities and network rearrangements at the microscopic level, by which the material deforms and flows. Such goals will be only partly achieved here, since, leaving the detailed study of velocity correlations and strain mechanisms for future work, we focus on simple characterizations that are local in space and time. We also check here that the various microstructural variables studied, if measured in D-type simulations, approach their values observed in S-type ones in the limit of (at least in the large system limit).

Packing geometry is classically described with a few state variables, among which the simplest ones are scalar: the volume fraction, the coordination number, as studied here in Section IV.1 below.

The much-studied distribution of contact force values Radjaï et al. (1996); Blair et al. (2001) is also determined in the present case (Section IV.2), and we check for effects of inertia and anisotropy.

Under stress, or influenced by the history of their assembling process, the microstructure of grain packings develops anisotropic features, which are most often characterized with the fabric tensors, expressing statistics on orientations of normal directions at contacts, as studied in Sec. IV.3. The critical state is microscopically characterized by stationary values of , , and fabric tensors, which are reached after a sufficiently large interval of monotonically growing strain in the quasistatic regime Radjaï et al. (2004); Radjaï and Roux (2004); Rothenburg and Kruyt (2004).

iv.1 Coordination number

The coordination number strongly depends on in steady state shear flow, and it is also affected by . It decreases with increasing , or with increasing . As to the influence of on , it is notable for the largest values explored, but it decreases as the quasistatic limit is approached. Larger viscous damping coefficients increase the duration of contacts in shear flow, and thus produce slightly better coordinated networks on average. However, the intensity of viscous forces becomes irrelevant in the quasistatic limit. According to our results, the dependence of can safely be ignored for . The dependence of is shown in Fig. 12.

Figure 12: (Color online) Coordination number as a function of inertial number , for (red square dots joined by dotted line, bottom data points) and (blue crosses, top, dashed line).

The coordination number of the equilibrated (S-type) anisotropic configurations is also very close to 6. This is a consequence of the isostaticity property of the force-carrying structure (also called backbone Agnolin and Roux (2007a)) of equilibrated sphere packings in the rigid limit – a remarkable property discussed in several recent publications O’Hern et al. (2003); Donev et al. (2005); Agnolin and Roux (2007a), which is specific to packings of rigid, frictionless and cohesionless spherical grains Roux (2000); Donev et al. (2007).

Fig. 12 shows that for quite low values of , many contacts are lost ( is down to about 5 for in the range). The proportion of rattlers increases with increasing : is less than in S simulations and for D simulations with and , but is equal to for and . Our results are compatible with the theoretical value in the limit and . Furthermore, it has been often observed that the grains only have a small number of contacts bearing large forces – this is the very reason why the “force chains” exist Ouaguenouni and Roux (1997); Radjaï et al. (1998); Somfai et al. (2005). Consequently, as contacts carrying smaller forces are necessarily shorter lived, and tend to rarefy as increases, the populations of grains with the largest local coordination are quickly depleted.

iv.2 Distribution of forces

Fig. 13 is a plot of the probability distribution function of the intergranular force normalized by the average force, for different values of . The force distribution strongly depends on : for the probability distribution function ( denoting the ratio of the normal force to the average value ) is monotonically decreasing. For smaller values of , has a maximum, around , and an approximately exponential decay for large values, as often observed in equilibrated granular packings Coppersmith et al. (1996); Radjaï et al. (1996); Ouaguenouni and Roux (1997); Blair et al. (2001); Silbert et al. (2002b); Donev et al. (2005); Agnolin and Roux (2007a).

Figure 13: (Color online) Probability distribution function of for (red triangles), (blue round dots) and (black square dots).

The distributions obtained for the low values of in D-type simulations gradually approach the one obtained in S-type, equilibrated packings under maximum shear stress. The Kolmogorov-Smirnov test Press et al. (1988) can be used to detect the influence of parameters on the force distribution – the answer depending of course on the level of statistics of the available data. Based on 10 independent configurations of 4000 grains, it leads to the conclusions that no significant difference in force distribution could be detected between S-type results under maximum shear stress and D-type ones, and no influence of either, provided the inertial parameter is small enough: , while some influence of is only visible for . Our results are also compatible with a unique distribution, valid for maximum shear stress equilibrium configurations as well as for isotropic ones.

iv.3 Fabric

Macroscopic friction is known to stem (at least partially) from the build-up of fabric anisotropy in materials made of frictional beads or disks Radjaï and Roux (2004); Radjaï et al. (2004). This connection is explored here with frictionless beads.

Anisotropy of the tridimensional contact network can be characterized by the probability density function of finding a contact with direction . is the colatitude angle and is the longitude angle of the spherical coordinates. can be expanded in a series of spherical harmonics. The coefficients of the expansion are in one-to-one correspondence with the values of the fabric tensors, which are defined as the moments of the distribution of normal unit vectors on the unit sphere. Since a contact is left invariant by the parity symmetry , satisfies . This means that the coefficients of odd order in the expansion in spherical harmonics are all equal to zero, and corresponds to the vanishing of all odd order fabric tensors. Coefficients can be computed from even order tensor products, viz.


denoting the set of contacts, labeled with index , where the normal unit vector is .

Keeping only the lowest order of anisotropy, the expansion of is restricted to the spherical harmonics of order two. Coefficients of the development are directly related to the value of the fabric tensor of order two, denoted by :


The constant corresponds to an isotropic distribution and the next five terms of the development characterize the anisotropy of the material at the lowest order. Functions are combinations of spherical harmonics of order two, with following expressions:


Fabric tensor is computed as a time average in the steady shear simulations. The numerical results show that and are always less than their respective statistical uncertainties, and can be considered as equal to zero, as requested by the symmetry in simple shear. We observe that is always greater (by at least one order of magnitude) than the two other non negligible anisotropic coefficients, and . These two latter terms are below for . Such low values are comparable with sample to sample fluctuations in equilibrated configurations. Thus, in the quasistatic limit, the anisotropy can be characterized by the sole coefficient, the limit of which, as , is evaluated as for and for , with a fitting procedure. Like and , strongly varies with (Fig. 14), and its dependence on can be represented by a power law, with exponent for and for .

Figure 14: (Color online) as a function of inertial number , with and , for (red square dots connected by a dotted line) and for (blue crosses connected by a dashed line). Both lines are power law fits of .

values measured in S-type equilibrated samples under maximum shear stress are influenced by system size , very similarly to the static friction coefficient: larger values of are observed (typically for ), but the excess over , the estimate from D-type simulations in the quasistatic limit, regresses as increases, and the extrapolated macroscopic limit is compatible with the estimated values of . This will be further examined in the more general context of the relationship between stress and anisotropies, for arbitrary stress tensors, in a forthcoming publication Peyneau and Roux (2008).

On changing from down to , increases (correlatively with the decrease of ), by about 30% for . This relative change is reduced to about 1% for and the effect of vanishes in the limit of .

Variations of with parameters , and are qualitatively understood on noting that is negatively correlated with the coordination number. If there are more contacting neighbors, on average, around a sphere, they are prevented by steric constraints from achieving highly anisotropic orientation distributions. This argument, which with simple assumptions was made quantitative in 2D in Ref. Radjaï et al. (2004), thus explains that the increase of observed as is lowered tends to reduce , as observed on Fig. 14. Similarly, the larger anisotropies observed away from the quasistatic limit are made possible by the smaller number of contacts. The increase of with is also due to the correlation of force intensities with contact directions: on evaluating separately the fabric of the subnetworks corresponding to forces larger (or smaller) than the average contact force, one typically obtains, for , values of twice as large (respectively: four times as small) as with the complete contact network. Contacts with small forces open if is increased, and the remaining more strongly loaded ones are consequently more anisotropically oriented.

V Discussion

This work was devoted to the study of frictionless identical spherical balls subjected to simple shear. The influence of the three dimensionless quantities controlling the problem – inertial number , stiffness number and level of viscous damping – was carefully assessed and we observed that has the most dramatic impact on the system behavior. Fluctuations of the measured quantities were shown to vanish for large systems. Consequently, the particular nature of the boundary conditions employed has no importance: for sufficiently large systems, fixed-volume simulations would lead to the very results we obtained with our stress driven numerical experiments. Particular attention was paid to the macroscopic geometric limit, that is the triple limit , and . In this régime, the system behavior is governed by a succession of instabilities due to dynamical rearrangements of the contact network. A thorough investigation of such events remains an interesting, yet challenging, perspective.

The existence of a nonzero macroscopic friction angle was evidenced by two different kinds of simulations – shear-rate controlled dynamic calculations (D-type simulations) and quasistatic stress-controlled calculations (S-type simulations). Whereas the dynamic friction angle is independent of the system size for , the static friction angle is very sensitive to the number of grains and is systematically greater than for all studied sizes () and increases for decreasing . This might be the reason why localization seems to occur more easily as the system size decreases. In finite-size systems, the shear stress is a multivaluated function of the strain rate in the quasistatic limit and the range of multivaluation increases with decreasing . Thus shear bands are more likely to appear in small systems da Cruz et al. (2005); Varnik et al. (2003, 2004). However, in the macroscopic geometric limit, we found that both friction angles and are equal within statistical uncertainties. In frictionless granular assemblies, all dissipation is due to viscous terms in contact forces, which therefore can be regarded Taboada et al. (2006) as the physical origin of macroscopic friction. However, the value of the damping coefficient is irrelevant in the quasistatic limit since the amount of dissipated energy is geometrically determined. In the macroscopic geometric limit, we have seen that the shear has no effect on the microscopic scalar quantities of the material (coordination number, distribution of forces), but it induces some structural anisotropy and a correlation between force intensities and contact orientations. We thus attribute the macroscopic friction angle to the shear-induced anisotropy of the material, as in the frictional case Azéma et al. (2007). Ref. Peyneau and Roux (2008) will show quantitatively that this is indeed the case.

The result that contrasts with observations on Lennard-Jones glasses at temperature  Varnik et al. (2003, 2004) and on granular avalanches Pouliquen (1999); Daerr and Douady (1999). Glass simulations show that the dynamic angle is less than the static one. This difference is linked to a stress overshoot visible on strain-stress curves. Similarly, in dense granular materials with friction, the shear stress goes through a maximum before the steady state (“critical state”) is reached, a feature which is absent in frictionless granular assemblies (both states coincide in this case). Similar differences () are reported for granular flows down inclined plane. Thus, in Refs. Pouliquen (1999); Daerr and Douady (1999), is less than , where is the inclination of the plane and the thickness of the flowing layer in the stationary state. The small thickness of the layer (typically less than ten grain diameters) and the intergranular friction are certainly responsible for this hysteresis.

The stress-dilatancy interplay is a well known feature of granular materials. However, our simulations show that homogeneously sheared frictionless bead assemblies do not display any dilatancy in the macroscopic geometric limit. In this limit, volume fraction remains equal to during the whole time the material is sheared and the backbone stays isostatic in the rigid and quasistatic limits. This surprising lack of dilatancy can be intuitively understood in the light of the simple model presented in Sec. III.5.3. We thus conclude that the steady state (critical) volume fraction is equal to .

The behavior of frictionless granular assemblies under arbitrary load directions will be the subject of a future work Peyneau and Roux (2008) in order to gain a better knowledge of the yield surface and of the mechanical properties of such granular systems under a small enough stress deviator (before failure).

One motivation of the present work is the study of highly concentrated non-Brownian suspensions (Péclet number Pe ), modeled as assemblies of nearly touching grains bonded by a viscous lubricant Melrose and Ball (1995); Ball and Melrose (1997); Melrose and Ball (2004a). Ideal lubrication effectively suppresses the tangent forces. Lubricated dynamics has already been employed as a means to obtain the force-carrying contact network of frictionless rigid particles, as the set of viscous bonds on which stresses concentrate Ouaguenouni and Roux (1995). Although crude, our current model should be able to reproduce the behavior of dense suspensions in the quasistatic limit. In this régime, the system evolves via a sequence of equilibrium states. At some point, the initial network is no longer able to sustain the imposed stress, it becomes unstable and a dynamic “crisis” occurs. Consequently, the evolution of the system is not quasistatic in the strictest sense (each point of the configuration space cannot be reached through a continuous series of equilibrium configurations). However, details of the dynamics are expected to be irrelevant. Thus we expect that the same equilibrium states will be visited in the quasistatic limit by both frictionless granular systems and dense suspensions with frictional grains. According to the simple toy model of Sec. III.5.3, a dense suspension might be sketched by a slider moving on a bumpy surface in a media of viscosity . Close to the quasistatic limit, the most important parameter would be the dimensionless number . One may notice that it is very similar to the parameter introduced by Cassar et al. that controls submarine avalanches in what they call the viscous regime Cassar et al. (2005). Steady shear simulations evidenced that the material is still able to flow with a volume fraction approximately equal to . This result is consistent with theoretical results pertaining to suspensions, where the volume fraction at which the viscosity of the suspension diverges is believed to tally with the random close packing volume fraction Sierou and Brady (2001). However, it is not in agreement with the experiments exposed in Ovarlez et al. (2006), where the value of was found to be below . This discrepancy very likely originates in small scale features of the experimental system that are not accounted for a model of perfectly lubricated spherical beads. The behavior of dense suspensions is known to be strongly impacted by short-range physics Melrose and Ball (2004b). In the near future, we plan to study lubricated pastes with frictional contacts in the spirit of the simplified Stokesian dynamics scheme proposed by Ball and Melrose Melrose and Ball (1995); Ball and Melrose (1997); Melrose and Ball (2004a).

Appendix A Crystallization under shear

Small samples, in both D and S-type simulations, tend to form strongly ordered structures under shear. This phenomenon, which do not occur for , is briefly reported here. A more detailed study would be outside the scope of the present paper, and would require some investigation of the role of cell shape and boundary conditions, which is necessarily important in such small systems.

2 out of 3 S-type samples with and stiffness level , and 2 out of 2 D-type samples with , and present the following anomalies. First, solid fractions are considerably higher than (and even more so considering the size effect O’Hern et al. (2003); Agnolin and Roux (2007a) on ), with values approximatively equal to (see fourth line of Tab. 3). Apparent friction coefficients are also particularly large. A lower bound for the static macroscopic friction coefficient of S-type ordered samples is , whereas dynamic macroscopic friction coefficient of D-type ordered samples for , , and may exceed by the corresponding friction coefficient in bigger samples that do not experience any ordering. S-type samples also have very large coordination numbers, . This latter characteristic is a clear indicator of partial crystalline order, as one necessarily has in generically disordered situations. The denser crystal arrangements, face-centered cubic (fcc) and hexagonal compact (hcp) (and stacking variants thereof), have . For D simulations, anomalous values of and appear after strains of order .

In order to detect crystalline order more quantitatively, we use the standard order parameters and employed in Volkov et al. (2002); Aste et al. (2005); Luchnikov et al. (2002); Agnolin and Roux (2007a). Values of the pair can be used to distinguish different local environments. In Agnolin and Roux (2007a), following Aste et al. (2005), the frequency of occurrence of ranges of values and , respectively corresponding to fcc-like and hcp-like configurations around one grain, were recorded. In the present case, most samples had very similar proportions of hcp-like and fcc-like local arrangements as in the RCP states studied in Agnolin and Roux (2007a): about 12% of beads fall in the hcp category, and fcc-like ones are virtually absent. The exceptions are the samples with anomalous, crystal-like properties, for which, while none of the beads has an fcc-like environment in that sense, the proportion of the hcp-like category raises to about 60% in S samples and to 40% in D ones.

A direct visualization, Fig. 15, reveals strikingly ordered configurations.

Figure 15: Crystalline order induced by shear in one S sample with .

A tentative conclusion to those preliminary observations is that the small samples tend to crystallize on somewhat shear-distorted hcp lattices. One convenient characterization of order that is not sensitive to the distortion of crystalline patterns was suggested in Volkov et al. (2002), and used in Agnolin and Roux (2007a). With this method, more than 90% of the particles of the anomalous samples are declared to belong to crystalline regions.


  1. Laboratoire des Matériaux et des Structures du Génie Civil 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.
  2. This numerical experiment displays shear-induced ordering, a feature observed only for systems of beads (see Appendix A for details).


  1. H. J. Herrmann, J.-P. Hovi, and S. Luding, eds., Physics of Dry Granular Media (Balkema, Dordrecht, 1998).
  2. H. Hinrichsen and D. E. Wolf, eds., The Physics of Granular Media (Wiley-VCH, Berlin, 2004).
  3. R. García Rojo, H. J. Herrmann, and S. McNamara, eds., Powders and Grains 2005 (Balkema, Leiden, 2005).
  4. A. Castellanos, Advances in Physics 54, 263 (2005).
  5. J. J. Stickel and R. L. Powell, Annu. Rev. Fluid Mech. 37, 129 (2005).
  6. G. Ovarlez, F. Bertrand, and S. Rodts, Journal of Rheology 50, 259 (2006).
  7. J. F. Brady, Journal of Chemical Physics 99, 567 (1993).
  8. D. Cumberland and R. Crawford, The Packing of Particles (Elsevier, Amsterdam, 1987).
  9. O. Reynolds, Philosophical Magazine (5th series) pp. 469–481 (1885).
  10. D. M. Wood, Soil Behaviour and Critical State Soil Mechanics (Cambridge University Press, 1990).
  11. P. W. Rowe, Proceedings of the Royal Society of London A269, 501 (1962).
  12. P. A. Cundall and O. D. L. Strack, Géotechnique 29, 47 (1979).
  13. J.-N. Roux and F. Chevoir, Bulletin des Laboratoires des Ponts et Chaussées 254, 109 (2005).
  14. J. F. Brady and G. Bossis, Annual Review of Fluid Mechanics 20, 111 (1988).
  15. F. da Cruz, S. Emam, M. Prochnow, J.-N. Roux, and F. Chevoir, Phys. Rev. E 72, 021309 (2005).
  16. T. Hatano, Physical Review E 75, 060301(R) (2007).
  17. A. Taboada, N. Estrada, and F. Radjaï, Physical Review Letters 97, 098302 (2006).
  18. C. Thornton, Géotechnique 50, 43 (2000).
  19. A. S. J. Suiker and N. A. Fleck, ASME Journal of Applied Mechanics 71, 350 (2004).
  20. D. Bideau and A. Hansen, eds., Disorder and Granular Media (Elsevier, 1993).
  21. C. O’Hern, L. E. Silbert, A. J. Liu, and S. R. Nagel, Physical Review E 68, 011306 (2003).
  22. A. Donev, S. Torquato, and F. H. Stillinger, Phys. Rev. E 71, 011105 (2005).
  23. I. Agnolin and J.-N. Roux, Phys. Rev. E 76, 061302 (2007a).
  24. G. Combe and J.-N. Roux, Physical Review Letters 85, 3628 (2000).
  25. J.-N. Roux and G. Combe, C. R. Académie des Sciences (Physique) 3, 131 (2002).
  26. A. V. Tkachenko and T. A. Witten, Physical Review E 62, 2510 (2000).
  27. J. D. Goddard and A. K. Didwania, Quaterly Journal of Mechanics and Applied Mathematics 51, 15 (1998).
  28. L. E. Silbert, D. Ertaş, G. S. Grest, T. C. Halsey, and D. Levine, Physical Review E 65, 031304 (2002a).
  29. N. Xu and C. S. O’Hern, Physical Review E 73, 061303 (2006).
  30. K. L. Johnson, Contact Mechanics (Cambridge University Press, 1985).
  31. E. Somfai, J.-N. Roux, J. Snoeijer, M. van Hecke, and W. van Saarloos, Phys. Rev. E 72, 021301 (2005).
  32. M. Allen and D. Tildesley, Computer simulations of liquids (Oxford University Press, Oxford, 1987).
  33. R. J. Bathurst and L. Rothenburg, Mechanics of Materials 9, 65 (1990).
  34. J. Christoffersen, M. M. Mehrabadi, and S. Nemat-Nasser, Journal of Applied Mechanics 48, 339 (1981).
  35. P. Rognon, J.-N. Roux, D. Wolf, M. Naaïm, and F. Chevoir, Europhysics Letters 74, 644 (2006).
  36. F. Da Cruz, F. Chevoir, D. Bonn, and P. Coussot, Phys. Rev. E 66, 051305 (2002).
  37. GDR MiDi, European Physical Journal E 14, 341 (2004).
  38. P. Jop, Y. Forterre, and O. Pouliquen, Nature 44, 727 (2006).
  39. I. Agnolin and J.-N. Roux, Phys. Rev. E 76, 061303 (2007b).
  40. E. Aharonov and D. Sparks, Phys. Rev. E 65, 051302 (2002).
  41. H. Flyvbjerg and H. Petersen, Journal of Chemical Physics 91, 461 (1989).
  42. F. Radjaï and S. Roux, in Hinrichsen and Wolf (2004), pp. 165–187.
  43. L. Rothenburg and N. P. Kruyt, International Journal of Solids and Structures 41, 5763 (2004).
  44. N. P. Kruyt and L. Rothenburg, Journal of Statistical Mechanics; Theory and Experiment p. P07021 (2006).
  45. L. Oger, S. Savage, D. Corriveau, and M. Sayed, Mechanics of Materials 27, 189 (1998).
  46. F. Radjaï, M. Jean, J.-J. Moreau, and S. Roux, Phys. Rev. Lett. 77, 274 (1996).
  47. D. L. Blair, N. W. Mueggenburg, A. H. Marshall, H. Jaeger, and S. R. Nagel, Physical Review E 63, 041304 (2001).
  48. F. Radjaï, H. Troadec, and S. Roux, in Granular Materials: Fundamentals and Applications, edited by S. J. Antony, W. Hoyle, and Y. Ding (Royal Society of Chemistry, Cambridge, 2004), pp. 157–183.
  49. J.-N. Roux, Physical Review E 61, 6802 (2000).
  50. A. Donev, R. Connelly, F. H. Stillinger, and S. Torquato, Phys. Rev. E 75, 051304 (2007).
  51. S. Ouaguenouni and J.-N. Roux, Europhysics Letters, 39, 117 (1997).
  52. F. Radjaï, D. E. Wolf, M. Jean, and J.-J. Moreau, Physical Review Letters 80, 61 (1998).
  53. S. N. Coppersmith, C. H. Liu, S. Majumdar, O. Narayan, and T. A. Witten, Physical Review E 53, 4673 (1996).
  54. L. E. Silbert, G. S. Grest, and J. W. Landry, Physical Review E 66, 061303 (2002b).
  55. W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, Numerical Recipes (Cambridge University Press, 1988).
  56. P.-E. Peyneau and J.-N. Roux (2008), in preparation.
  57. F. Varnik, L. Bocquet, J.-L. Barrat, and L. Berthier, Phys. Rev. Lett. 90 (2003).
  58. F. Varnik, L. Bocquet, and J.-L. Barrat, J. Chem. Phys. 120 (2004).
  59. E. Azéma, F. Radjaï, R. Peyroux, and G. Saussine, Phys. Rev. E 76 (2007).
  60. O. Pouliquen, Phys. Fluids 11, 542 (1999).
  61. A. Daerr and S. Douady, Nature 399, 241 (1999).
  62. J. Melrose and R. Ball, 32, 535 (1995).
  63. J. Melrose and R. Ball, J. Rheol. 48, 961 (2004a).
  64. R. Ball and J. Melrose, Physica A 247, 444 (1997).
  65. S. Ouaguenouni and J.-N. Roux, Europhys. Lett. 32, 449 (1995).
  66. C. Cassar, M. Nicolas, and O. Pouliquen, Phys. Fluids 17 (2005).
  67. A. Sierou and J. Brady, J. Fluid Mech. 448, 115 (2001).
  68. J. Melrose and R. Ball, J. Rheol. 48, 937 (2004b).
  69. T. Aste, M. Saadatfar, and T. J. Senden, Phys. Rev. E 71 (2005).
  70. I. Volkov, M. Cieplak, J. Koplik, and J. R. Banavar, Physical Review E 66, 061401 (2002).
  71. V. Luchnikov, A. Gervois, P. Richard, L. Oger, and J.-P. Troadec, Journal of Molecular Liquids 96-97, 185 (2002).