A Rigidity and stiffness matrices

Computer simulation of model cohesive powders: influence of assembling procedure and contact laws on low consolidation states.


Molecular dynamics simulations are used to investigate the structure and mechanical properties of a simple two-dimensional model of a cohesive granular material. Intergranular forces involve elasticity, Coulomb friction and a short range attraction akin to the van der Waals force in powders. The effects of rolling resistance (RR) at intergranular contacts are also studied. The microstructure of the cohesive packing under low pressure is shown to depend sensitively on the assembling procedure which is applied to the initially isolated particles of a granular gas. While a direct compression produces a final equilibrated configuration with a similar density to that of cohesionless systems, the formation of large aggregates prior to the application of an external pressure results in much looser stable packings. A crucial state variable is the ratio of applied pressure , acting on grains of diameter , to maximum tensile contact force . At low the force-carrying structure and force distribution are sensitive to the level of velocity fluctuations in the early stages of cluster aggregation. The coordination number of packings with RR approaches 2 in the limit of low initial velocities or large rolling friction. In general the force network is composed of hyperstatic clusters, typically comprising 4 to a few tens of grains, in which forces reach values of the order of , joined by barely rigid arms, where contact forces are very small. Under growing , it quickly rearranges into force chain-like patterns that are more familiar in dense systems. Density correlations are interpreted in terms of a fractal structure, up to a characteristic correlation length of the order of ten particle diameters for the studied solid fractions. The fractal dimension in systems with RR coincides, within measurement uncertainties, with the ballistic aggregation result, in spite of a possibly different connectivity, but is apparently higher without RR. Possible effects of micromechanical and assembling process parameters on mechanical strength of packings are evoked.


81.05.Rm: Porous materials, granular materials,
83.10.Rs: Computer simulation of molecular and particle dynamics,
61.43.Hv: Fractals, macroscopic aggregates,
47.57.J-: Colloidal systems

I Introduction

Granular materials are currently being studied by many research groups Herrmann et al. (1998); Kishino (2001); Hinrichsen and Wolf (2004); García Rojo et al. (2005), motivated by fundamental issues (such as the relations between microstructure and global properties) as well as by practical needs in civil engineering and in the food and drug industries. The relation of their mechanical behavior in quasistatic conditions to the packing geometry, which depends itself on the assembling procedure, tends to escape intuition and familiar modelling schemes.

The configuration of the contact networks is hardly accessible to experiments, even though particle positions are sometimes measured Richard et al. (2003); Xu et al. (2004); Aste et al. (2004, 2005), and some experimental quantitative studies on intergranular contacts carried out in favorable cases (such as millimeter-sized beads joined by capillary menisci Xu et al. (2004); Kohonen et al. (2004); Fournier et al. (2005); Richefeu et al. (2006)). Intergranular forces are also, most often, inaccessible to measurements. Consequently, computer simulation methods of the “discrete element” type, as introduced 30 years ago Cundall and Strack (1979), have proved a valuable tool to investigate the internal states of granular systems. Simulation methods like molecular dynamics Herrmann and Luding (1998) or “contact dynamics” Jean (1999); Jean et al. (2001) have been gaining an increasingly large constituency of users and a wide range of applications, as witnessed e.g., by recent conference proceedings García Rojo et al. (2005).

Dry assemblies of grains interacting via contact elasticity and friction, such as sands or glass beads, might form stable packings of varying solid fraction (typically between 58% and 64% for monosized spheres if they do not crystallize), which deform plastically in response to changes in stress direction, rather than stress intensity. Their elastic or elastoplastic properties have been studied by discrete simulation (see, e.g., Thornton (2000); Suiker and Fleck (2004)), and, in agreement with laboratory experiments and macroscopic modelling Wood (1990), found to depend sensitively on the initial density. Numerical simulation also stressed the importance of additional variables such as coordination number Agnolin and Roux (2005) and fabric Bathurst and Rothenburg (1990); Radjai and Roux (2004), and it has often been applied to the study of quasi-static stress-strain behavior of granular assemblies (refs. Thornton (2000); Roux and Combe (2002); Radjai and Roux (2004); Suiker and Fleck (2004) are a few examples among a large literature).

Cohesive grains exhibit much larger variations in their equilibrium densities, and they are sensitive to stress intensity as well as direction : on increasing the confining pressure, the specific volume of a clay can irreversibly decrease by a factor of 4 Mitchell (1993). Likewise, series of experiments carried out in the Seville group on model powders Castellanos (2005) (xerographic materials) in which the strength of van der Waals attraction is controlled by additives covering part of the grain surfaces, reveal a similar variation of porosity with confining pressure. It is notable that such packings of particles of rotund shape and nearly the same size can stay in mechanical equilibrium at much lower solid fractions (down to 25-30%) than cohesionless granular systems.

Despite this wider variety of equilibrium structures and mechanical behaviors, cohesive granular materials have much less frequently been investigated by numerical simulation than cohesionless ones.

Some of the recent numerical studies, such as those of refs. Yang et al. (2000, 2003a) have investigated the packing structures of spherical beads deposited under gravity, depending on micromechanical parameters, including adhesion strength. Another set of publications report on simulations of the dynamical collapse and compaction, both in two Kadau et al. (2002, 2003); Wolf et al. (2005) and three Bartels et al. (2005) dimensions, the main results being the relations between density and pressure increments, and their dependence on micromechanical parameters. Some works focussed on the fracture of bound particle assemblies in static Delenne et al. (2002, 2004) or dynamic Thornton and Liu (2004) conditions, others on wet bead packs in which cohesion stems from liquid bridges joining neighboring particles, investigating the structure of poured samples Yang et al. (2003b) or the shear strength Richefeu et al. (2006) of such materials. These two latter types of studies deal with relatively dense materials, as does the numerical biaxial compression test of Luding (2005). Flow of cohesive materials has also been addressed in recent publications Rognon et al. (2005, 2006); Brewster et al. (2005).

Yet, numerical studies of the mechanics of loose, solid-like cohesive granulates are quite scarce. This contrasts with the abundant literature on the geometry of model loose particle packings and colloidal aggregates, which tend to form fractal structures. Refs. Smirnov (1990); Meakin (1999) are useful overviews of aggregation processes and the geometric properties of the resulting clusters, as obtained by numerical simulation. In such processes, particle aggregates are usually regarded as irreversibly bound, rigid solids, while the interaction between separate clusters reduces to a “sticking rule”, so that both intra- and inter-aggregate mechanical modelling is bypassed. Interestingly, one simulation study Bratberg et al. (2002) shows that structures resulting from geometric deposition algorithms are not always stable once a mechanical model is introduced.

It seems that numerical simulations of both geometric and mechanical properties of loose granular assemblies forming solid aggregates are still lacking.

The present paper addresses part of this issue. It reports on numerical simulation studies of cohesive granular materials, with the following specificities:

  • the assembling process is simulated with the same mechanical model as applied to solid-like configurations, and its influence on the packing microstructure is assessed ;

  • special attention is paid to loose particle packings in equilibrium under vanishing or low applied pressure;

  • both geometric and mechanical properties are investigated ;

  • isotropic and homogeneous systems are studied, as representative samples for bulk material properties.

We consider a simple model system in two dimensions, introduced in section II, along with the numerical simulation procedure. Despite its simplicity we shall see that this model yields results that are amenable to comparisons with experimental situations.

Section. III is devoted to the important issue of the procedure to prepare samples, and its influence, as well as that of micromechanical features such as rolling resistance (RR), on final density and coordination number in solid packings in equilibrium. In Section IV we investigate the force distributions and force patterns of the equilibrated loose configurations under vanishing or low applied pressure. Some specific aspects of the force-carrying structures in low density assemblies will be studied and related to the assembling process. In Section V, we characterize the geometry and density correlations in loose samples, resorting to the fractal model traditionally employed for colloidal aggregates. Finally we conclude in Section VI with a few remarks about future improvements and further developments of this work, some of which will be presented in a forthcoming publication Gilabert et al. ().

Ii Model material

ii.1 System definition, equations of motion

We consider a two-dimensional model material: an assembly of disks with diameters uniformly distributed between and . The maximum diameter, , will be used as unit of length. The mass of grain is and its moment of inertia , i.e. disks are regarded as homogeneous bodies and the mass of a disk of maximum diameter is the unit of mass.

The disks are enclosed in a rectangular cell the edges of which are parallel to the axes of coordinates and , with respective lengths and . Periodic boundary conditions are used, thereby avoiding wall effects. Neighboring grains, say and , might interact if they are brought into contact or very close to each other, hence a force and a moment exerted by onto at the contact point. Simulations do not model material deformation in a contact region, but consider overlapping particles, and the contact point is defined as the center of the intersecting surface of the two disks. In the case of an interaction without contact, the force will be normal to the surfaces at the points of nearest approach, and therefore carried by the line of centers. Let denote the position of the center of disk . is the vector joining the centers of and , and their overlap distance. The degrees of freedom, in addition to the positions , are the angles of rotation , velocities , angular velocities of the grains (), the dimensions of the cell containing the grains and their time derivatives, through the strain rates:

in which denotes the initial size for the corresponding compression process. The time evolution of those degrees of freedom is governed by the following equations.


In Eqns. (1) and (2), only those disks interacting with , i. e. in contact or very close, will contribute to the sums on the right-hand side. In Eqn. (3), is the externally imposed stress component, is the measured stress component, resulting from ballistic momentum transport and from the set of intergranular forces , denotes the cell surface area, and is a generalized inertia parameter.

Stresses and , rather than strains or cell dimensions, are controlled in our simulation procedure. Note that compressions are counted positively for both stresses and strains. Eqn. (3) entails that the sample will expand (respectively, shrink) along direction if the corresponding stress is larger (resp., smaller) than the requested value , which should be reached once the system equilibrates. This barostatic method is adapted from the ones initially proposed by Parrinello and Rahman Parrinello and Rahman (1980, 1981, 1982) for Hamiltonian, molecular systems.

The choice of the “generalized mass” is rather arbitrary, yet innocuous provided calculations are restricted to small strain rates. In practice we strive to approach mechanical equilibrium states with good accuracy, and choose in order to achieve this goal within affordable computation times. We usually attribute to a value equal to a fraction of the sum of grain masses (3/10 in most calculations), divided by a linear size of the cell. This choice is dimensionally correct and corresponds to the appropriate time scale for strain fluctuations in the case of a thermodynamic system.

ii.2 Interaction law

The contact law in a granular material is the relationship between the relative motion of two contacting grains and the contact force. As we deal with particles that may attract one another at short distance without touching, the law governing intergranular forces and moments is best referred to simply as the interaction law.

Although the interaction we adopted is based on the classical linear “spring-dashpot” model with Coulomb friction for contact elasticity, viscous dissipation and sliding, as used in many discrete simulations of granular media Herrmann and Luding (1998); Silbert et al. (2002); Brewster et al. (2005); Rognon et al. (2006), some of its features (short-range attraction and rolling resistance) are less common ; moreover, one can think of different implementations of the Coulomb condition, depending on which parts of normal and tangential force components are taken into account. Therefore, for the sake of clarity and completeness, we give a full, self-contained presentation of the interaction law below.

We express intergranular forces in a mobile system of coordinates with axes oriented along the normal unit vector (along ) and the tangential unit vector ( is a direct base in the plane), and use the convention that repulsive forces are positive.

The intergranular force , exerted by grain onto its neighbor is split into its normal and tangential components, thus defining scalars and . comprises a static term depending on the distance between disk centers, combining contact elasticity and distant, van der Waals type attraction, as shown on Fig. 1a, and a velocity-dependent viscous term . (Fig. 1b) is due to the tangential elasticity in the contact, and is limited by the Coulomb condition.

Figure 1: Graphical representation of the model for the adhesive elastic contact force as a function of the distance between the surfaces of particles and , . (a) The elastic normal force consists of a repulsive Hookean part plus a linearized attractive part . (b) The elastic tangential force is limited by the Coulomb cone (adhesion shifting its tip to on the normal force axis).

If disks and are not in contact, both the tangential component of force and the viscous part of the normal component vanish, while and still attract each other if the gap () between their surfaces is smaller than the attraction range ():


This expression is a linear approximation of a realistic van der Waals force law (see Fig. 1a), and contains two essential parameters, maximum attractive force , range, . Typically, is of the order of , being a superficial energy, the typical size of asperities Forsyth and Rhodes (2000) and is in the nanometer range.

In the case of contacting disks (), the attractive term is kept constant, equal to , while strains in the contact region result in normal () and tangential () elastic forces. It is also assumed that a viscous normal term opposes relative normal displacements. One thus writes:


The different terms introduced in Eqn. (5) are defined according to the following models. First,

is the linear elastic unilateral repulsion, due to the normal deflection in the contact as the disks are pressed against each other. is the normal stiffness coefficient, related to the elastic moduli of the material the grains are made of.

The viscous normal force opposes the normal relative receding velocity as long as the contact persists. The relative normal motion of two disks and in contact is that of an oscillator with viscous damping, and is the damping coefficient. We choose its value as a constant fraction of the critical damping coefficient,


This is equivalent to the choice of a constant restitution coefficient in normal collisions if . In the presence of attractive forces the apparent restitution coefficient in a collision will depend on the initial relative velocity, and will be equal to zero for small values, when the receding velocity after the collision will not be able to overcome the attraction and separate the particles. The minimum receding velocity for two particles of unit mass (i.e., of maximum diameter ) to separate is , with


The elastic tangential force in contact is linearly related to the elastic part of the total relative tangential displacement , as

and is subject to the Coulomb inequality. is the tangential stiffness coefficient. can be updated for all closed contacts according to

and vanishes as soon as the contact opens. Its elastic part satisfies

in which denotes the Heaviside function. This last equation introduces the friction coefficient . It is important to note that the Coulomb inequality,


applies to the sole repulsive elastic component of the normal force (see Fig. 1b). We chose not to implement any tangential viscous force.

The moment that disk exerts onto its contacting neighbor , of radius , in its center, is denoted by in Eqn. (2). It is first due to the tangential contact force, then to a possible moment of the force density distribution within the contact region. One thus writes:


is most often neglected on dealing with smooth, convex particle shapes, because the contact region is very small on the scale of the particle radius.

To model rolling resistance (RR), like in Tordesillas and Stuart (2002), we introduce a rotational stiffness parameter and a rotational friction parameter in contacts, so that rolling elasticity and rolling friction are modelled just like sliding elasticity and friction. One thus writes

while enforcing the inequality


This involves the definition of as the elastic part of the total relative rotation . The total relative rotation angle satisfies

while the equation for is

Parameters and are often related to the size of a contact region Kadau et al. (2002). is dimensionally the product of a stiffness by the square of a length, which is of the order of the contact size. In the following we set to , while , which has the dimension of a length, is chosen equal to . The motivation for the introduction of RR into our model is twofold. First, cohesive particles are usually small (typically less than m in size) and irregular in shape. Contacts between grains are likely to involve several asperities, and hence some lateral extension, of the order of the distance between asperities, however small the normal deflection . Then, it will be observed that even quite a small rotational friction has a notable influence on the microstructure of cohesive packings.

ii.3 Control parameters and dimensional analysis

In this section we present the dimensionless parameters which express the relative importance of different physical phenomena. Such parameters enable qualitative comparisons with real materials, bearing in mind that the present model is admittedly an idealization of real powders and that our simulations do not aim at quantitative accuracy.

Dimensionless numbers related to contact behavior are the reduced interaction range , the friction coefficient , the viscous damping parameter , and the stiffness parameter .

Under the attractive force , the elastic deflection of one contact is


The stiffness parameter characterizes the amount of elastic deflection under contact force , relative to grain size (). A suitable analogous definition for Hertzian spheres in three dimensions would be .

The dimensionless number is the ratio of elastic to adhesive stiffnesses, and its physical meaning is similar to that of the Tabor parameter  Maugis (2000) for a Hertzian contact between spheres of diameter when the material Young modulus is and the interfacial energy is (more precisely, the equilibrium normal deflection , due to adhesion, in the contact between an isolated pair of grains, satisfies in this case).

The viscous damping parameter, , corresponds to a normal restitution coefficient in the absence of cohesion ().

In our calculations we set , corresponding to a high viscous dissipation in collisions, or a very low restitution coefficient in binary collisions. Models with a constant were adopted in other published simulation works Silbert et al. (2002); Somfai et al. (2005), although little is known about dissipation in collisions. is known to influence the packing structures obtained in the initial assembling stage Zhang et al. (2001); Silbert et al. (2002), but we did not investigate its effects in the present study. The simulations reported in Yang et al. (2000, 2003a) use the viscous force model introduced in ref. Brilliantov et al. (1996), with a choice of parameters corresponding to strongly overdamped dynamics (i.e., analogous to in our case).

In addition to those control parameters determined by the contact behavior, other dimensionless numbers are introduced by the loading or the process being applied to the material. The effect of the external pressure, compared to the adhesion strength, is characterized by a dimensionless reduced pressure :


In the present paper, we focus on the assembling process and the low range. As we shall see below (Sec. III) low density, tenuous structures are then stabilized by adhesion, and the relevant force scale is . However, as briefly reported in Gilabert et al. (2005), such structures tend to collapse upon increasing . These phenomena will be the subject of another paper Gilabert et al. (). Wolf et al. Wolf et al. (2005) introduced a dimensionless stress proportional to , and observed, in numerical simulations, stepwise increases in pressure to produce large dynamical collapse effects around . The importance of was also stressed in simulations of cohesive granular flow, in which the effects of cohesion on rheological laws were expressed in terms of a cohesion number defined as  Rognon et al. (2006). In three dimensions, should be defined as .

For large reduced pressures, externally imposed forces dominate the adhesion strength, and one should observe behaviors similar to those of confined cohesionless granular materials. For , the relevant force scale is . The influence of , which should then be defined as , so that the typical contact deflection satisfies , was studied in simulations of grains without adhesion Combe and Roux (2003). Whatever the reference force used to define it, the limit of rigid grains is . With relatively soft grains (say, below ), a significant number of additional contacts appear in dense configurations, due to the closing of gaps between near neighbors. Such a parameter defined with reference to pressure, in the case of contacts ruled by Hertz’s law between spherical grains made of a material with Young modulus , should be chosen as , in order to maintain .

In order to stay within the limit of rigid grains both for small and large , we choose quite a large value of : or .

Table 1 summarizes the values (or the range of values) of dimensionless parameters in the simulations presented below.

, , 0, ,
Table 1: Values of dimensionless model parameters used in most simulations. Note that is fixed by and to or . In the absence of cohesion, or for values of , is defined as

In addition to those values of the parameters, adopted as a plausible choice for realistic orders of magnitudes, some calculations were also performed with deliberately extreme choices, such as very large RR () or absence of friction ( and ), in order to better explore some connections between micromechanics and macroscopic properties. The corresponding results will be described in Section IV.

The definition of dimensionless parameters, suitably generalized to three-dimensional situations as and for spherical particles of diameter , enables one to discuss qualitative features and orders of magnitude in the model system defined with the parameters of Table 1 with comparisons to some cohesive packings studied in the laboratory.

When adhesive forces are due to liquid menisci joining neighboring particles, we should take , where is the surface tension. corresponds then to confining pressure in the range of 10-100 Pa for millimeter-sized particles, taking standard values for . Those are rather low pressures in practice, which are comparable, e.g., to the ones caused by the weight of a typical laboratory sand sample. Thus wet granular materials are commonly under reduced pressures of order 1 or larger, and are not observed with much lower solid fractions than dry ones Xu et al. (2004); Kohonen et al. (2004); Fournier et al. (2005); Richefeu et al. (2006).

The cohesive powders studied in refs. Watson et al. (2001); Sánchez Quintanilla (2003); Castellanos et al. (2005); Castellanos (2005) are xerographic toners with typical particle diameter m. , the van der Waals attractive force, is a few tens of nN, and the range is several nanometers Krupp (1967). Therefore, a reduced pressure would correspond to about  Pa in the experimental situation Valverde et al. (2004). This is an initial state of very low consolidation stress, which is present in a powder under gravity, provided a controlled gas flow, going upwards through the powder, counterbalances part of its weight Valverde et al. (2004). As to contact stiffnesses, our values of would correspond to  GPa (for ) or  GPa (for ), while the ratio would imply an interaction range of 10 nm. This gives us the correct orders of magnitudes for the toner particles, those being made of a relatively soft solid (polymer, such as polystyrene) with  GPa. Xerographic toner particles appear to undergo plastic deformation in the contacts Sánchez Quintanilla (2003); Maugis and Pollock (1984); Quintanilla et al. (2001); Gilabert et al. (2006). Plastic deflections of contacts are accounted for in the model of ref. Luding (2005), applied to the simulation of a biaxial compression of a dense powder. In our study, for simplicity’s sake, and because we expect macroscopic plasticity of loose samples to be essentially related to the collapse of tenuous structures, we ignored this feature.

ii.4 Equilibrated states

Although numerical simulations of the quasistatic response of granular materials requires by definition that configurations of mechanical equilibrium should be reached, equilibrium criteria are sometimes left unspecified, or quite vaguely stated in the literature. Yet, in order to report results on important, often studied quantities like the coordination number or the force distribution, it is essential to know which pairs of grains are in contact and which are not. Due to the frequent occurrence of small contact force values, this requires forces to balance with sufficient accuracy. We found that the following criteria allowed us to identify the force-carrying structure clearly enough. We use the typical intergranular force value to set the tolerance levels. A configuration is deemed equilibrated when the following conditions are fulfilled:

  • the net force on each disk is less than , and the total moment is lower than ;

  • the difference between imposed and measured pressure is less than ;

  • the kinetic energy per grain is less than .

We observed that once samples were equilibrated according to those criteria, then the Coulomb criterion (8), as well as the rolling friction condition (10) were satisfied as strict inequalities in all contacts. No contact is ready to yield in sliding, and with RR no contact is ready to yield in rolling either.

Iii Assembling procedure

It has been noted in experiments Mitchell (1993) and simulations Thornton (2000); Agnolin and Roux (2005); Silbert et al. (2002) that the internal structure and resulting behavior of solid-like granular materials is sensitive to the sample preparation procedure, even in the cohesionless case.

In the case of powders, it has been observed that the sedimentation in dry nitrogen (to minimize the capillary effects of the humidity on the interparticle adhesion) of a previously fluidized bed produced reproducible states of low solid fractions (down to ) Valverde et al. (2000); Castellanos et al. (2001). This initial state under such a low consolidation, as we commented in II.3, plays a decisive role on the evolution of the dynamics of powder packing. That is, appreciable differences in initial states will lead to considerable ones in final packings Castellanos et al. (2005). This is mainly due to the role of aggregation, which we shall analyze in the second part of this section.

The motivation of this section is to investigate the dependence on packing procedure in a cohesive granular system, the first step being to obtain stable equilibrated configurations with low densities. For comparison, some simulation results are presented for the same model material with no cohesion.

Specimens were prepared in two different ways, respectively denoted as method 1 and method 2, and the resulting states are classified as type 1 or type 2 configurations accordingly.

Due to our choice of boundary conditions, our samples will be completely homogeneous, under a uniform (isotropic) state of stress. This choice is justified by the complexity of seemingly more “realistic” processes, such as gravity deposition, due to the influence of many material (such as viscous dissipation, as recalled in Section II.3) and process parameters. Both pouring rate and height of free fall should be kept constant during such a pluviation process in order to obtain a homogeneous packing Zhang et al. (2001); Emam et al. (2005) with cohesionless grains. Cohesive ones, because of the irreversible compaction they undergo on increasing the pressure, would end up with a density increasing with depth. Hence our choice to ignore gravity in our simulations. Our final configurations should be regarded as representative of the local state of a larger system, corresponding to a local value of the confining stress.

iii.1 Method 1

In simulations of cohesionless granular materials, a common procedure Thornton (2000); Makse et al. (2004); Suiker and Fleck (2004) to prepare solid samples consists in compressing an initially loose configuration (a “granular gas”), without intergranular contacts, until a state of mechanical equilibrium is reached in which interparticle forces balance the external pressure (further compaction being prevented by the jamming of the particle assembly). We first adopted this traditional method, hereafter referred to as method 1, to assemble cohesive particles.

In this procedure, disks are initially placed in random non-overlapping positions in the cell, with zero velocity. We denote such an initial situation as the -state. Then the external pressure is applied, causing the cell to shrink homogeneously. Thus contacts gradually appear and the configuration rearranges until the system equilibrates at a higher density.

Examples of equilibrated configurations are shown on Fig. b, with and without cohesion.

(a) No cohesion,
(b) Cohesive system, ,
Figure 2: (Color online) Aspect of force-carrying structures in cohesionless and cohesive samples. Contact forces are displayed with the usual convention that the width of the lines joining the centers of interacting pairs of disks is proportional to the normal force, on scale (left) and (right). Red, green, and blue lines distinguish compressive, tensile, and distant interactions in the cohesive case, while rattlers appear in grey in the cohesionless sample.

This state is characterized by its solid fraction () and its coordination number , defined as the average number of interactions (contacts and distant attractions) for a particle in the packing, when the applied pressure is significantly smaller than (), as in the case of small powder samples assembled under gravity. With the values indicated above (at the end of Section II.3) for toner particles, (the relevant pressure scale in 3D) is of the order of  Pa, which corresponds to a normal consolidation stress in a cohesive powder with solid fraction Valverde et al. (2004).

In the absence of cohesion, the value of the applied pressure does not affect the properties of the packing (apart from setting the scale of intergranular forces) provided the typical contact deflection, , is small enough (rigid particle limit). We set this ratio to the value of in the cohesive case, i.e., equal to (see table 1), so that typical contact forces are of the same order of magnitude (due either to or predominantly to ) in both cases.

Effects of the initial solid fraction in the -state, and of cohesion, friction and rolling resistance parameters on and were measured in 3 sets of samples, with , 0.36 and . Each set consisted of configurations with the initial disorder (particle radii and initial positions) abiding by the same probability distribution, and the same number of particles (). The values of the friction coefficient, , used in these tests were 0.15 and 0.5. The values of and in these samples are listed in tables 2 and 3. Each one is an average on the different samples, and the indicated uncertainly is equal to the standard deviation.

Non-Cohesive samples
no RR RR
Solid fraction
Coordination number
Table 2: Solid fractions and coordination numbers obtained at the preparation of the specimens in equilibrated samples under for non-cohesive particles, using method 1.
Cohesive samples
no RR RR
Solid fraction
Coordination number
Table 3: Solid fractions and coordination numbers obtained at the preparation of the specimens in equilibrated samples under , for cohesive particles, using method 1.

Tables 2 and 3 show that the introduction of cohesion reduces the solid fraction at equilibrium, but this is a limited effect (less than 10% density reduction), which is quite insufficient to account for experimental observations. Unlike powders or clays, 2D particle packings with cannot undergo very large plastic density increases.

Theses tables also show that the increase of the friction coefficient and/or the inclusion of rolling resistance in the model tend to hinder motions and stabilize looser, less coordinated configurations, which results in a decrease of and .

However, the observed differences are rather small, especially in cohesionless systems. To evaluate the influence of RR with a given value of , we define , as the relative average decrease of the quantity due to the existence of RR. For example, for results differ by a mere and , and for , these variations are and . Comparing the effect of RR on with for and , one has . Likewise, for coordination numbers , one observes . This shows a clear correlation of variations introduced by friction and RR. The data in the non-cohesive case also exhibit very little dependence on initial density .

Results on cohesive systems show similar variations with the parameters of the contact model (friction and RR), but depend somewhat more sensitively on .

More refined information on the contact network is provided by the distribution of local coordination numbers, i.e. the proportions of particles interacting with neighbors, for (higher values were not observed).

Figure 3: Distribution of local coordination numbers (percentage of total particle number), without (left graph) and with (right graph) cohesion.

This distribution is depicted on Figure 3, for both cohesive and non-cohesive samples. These results gather information from all the statistically equivalent simulated samples, and slight corrections were applied in order to ignore the contacts with “rattler” particles in the non-cohesive case. Such particles are those that are free to move within the cage of their near neighbors, and transmit no force once the system is equilibrated. If they happen to be in contact with the backbone (i.e., the force-carrying structure), then the forces carried by such contacts should be below the tolerance set on the equilibrium requirement, and can safely be ignored. This is how the population of rattlers is identified. We observe that it can involve up to 18% of the total number of grains in the absence of cohesion (see Fig a).

This contrasts with the cohesive case, for which nearly all the grains are captured by the force-carrying structure because of attractive forces, and the rattlers are virtually absent. The particles with one contact equilibrate when the deflection of that contact is , as defined in (11). With RR, such a particle is entirely fixed. Without RR, it is only free to roll without sliding on the perimeter of its interacting partner, because such a contact is able to transmit a tangential force smaller than or equal to .

Without cohesion, the coordination of the force-carrying structure can be characterized with a coordination number , different from :


is the average number of contacts bearing non-negligible forces per particle on the backbone. Without cohesion, the backbone (or set of non-rattler grains) is the rigid part of the packing. With cohesion and RR, the whole interacting contact network is to be considered in order to study the rigidity properties of the system, and there are nearly no particles to eliminate. With cohesion and no RR, we observe in the samples obtained by the presently employed procedure (method 1) that the network of interparticle contacts or interactions is also rigid, apart from the free rolling of isolated grains with only one contact. (The rigidity properties of equilibrated samples are discussed below in Section IV and Appendix A).

Cohesive samples in equilibrium also comprise a small number of pairs of particles interacting without contact, i.e. separated by a gap smaller than the range of attraction, . These are only a small fraction, below 1%, of interacting pairs. Such pairs do not contribute to dissipation, since the frictional and viscous force components are only present in true contacts between neighboring grains. We observed that the time necessary to equilibrate the sample tend to increase when such distant interacting pairs are more numerous.

In addition to the elimination of free rattlers, the most notable effect of cohesion on local coordination numbers (Fig. 3) is to increase the proportion of disks with 2 contacts. Without cohesion, the Coulomb condition restricts the angle between the directions of the 2 contacts to values between and , where is the friction angle (). Thus if is small, a disk with two contacts should have its center close to the line of centers of its two partners. The increase of the population of 2-coordinated disks as is raised from to (see Fig. 3) in cohesionless systems corresponds to a less severe geometric restriction on contact angles. With cohesion, contacts may transmit a tangential force reaching while the normal force is equal to zero. Consequently, a disk might be in equilibrium with two contact points in arbitrary positions on its perimeter. As there is no geometric constraint on the angle between the two contact directions, 2-coordinated disks are easier to stabilize, and their proportion raises from about without cohesion to above with cohesion in the case . A population of disks with one contact (therefore carrying a vanishing normal force, with deflection ) is also present. Those particles are fixed by a small rolling resistance, but are free to roll on their interacting neighbor without RR. Such a rolling motion is not damped in our model. Therefore, on waiting long enough, they should eventually stop after a collision, in a stable position with 2 contacts. Such a collision is bound to happen because the contact network is completely connected. However, we stop our calculations when the kinetic energy is below a set tolerance (see Section II.4), and we do not wait until all freely rolling disks reach their final position. Hence the remaining population of disks with one contact in samples without RR.

The final configuration, with this preparation method, depends somewhat on the rate of compaction in the assembling stage. The latter is related to the choice of the dynamical parameter , the “mass” with which the changes in cell dimensions are computed with Eqn. (3). The slight influence of the initial solid fraction, , also relates to such dynamical effects: a lower value of entails larger colliding velocities, which favors larger final solid fractions.

Although some of the aspects of the model (in particular the homogeneous shrinking imposed through the periodic cell dimensions in a dynamical regime) do not correspond to experimental conditions, configurations of type 1 should be regarded as typical results of fast assembling processes, in which the particles are requested to balance the external pressure before stable loose structures can be built. When the toner particles mentioned at the end of Section II.3 are first fluidized, and then settle under their own weight, a rough estimate of the settling time, assuming particles are settling individually in air, and fall over distances of order 1 cm, is  s. Fig. 4, with the value  s corresponding to such particles, shows that the duration of the “method 1” compression process is a few milliseconds. In practice, due to the presence of the surrounding fluid, the packing of a powder in a loose state by settling and compaction of an initially fluidized state is therefore considerably slower than this numerical process.

In the next section, we consequently turn to the opposite limit, in which the external confining pressure is felt only after large, tenuous contact networks are formed.

iii.2 Method 2

Numerical procedure.

Figure 4: Solid fraction versus time for both preparation procedures, showing some aspects of the configurations at different stages. Point A is the initial state (or ). Aspects of configurations are shown for intermediate states B and B, and for final equilibrated states C and D (at ). Point C corresponds to the stage when all disks are assembled in a unique aggregate, then equilibrated at (both aggregation and equilibration stages take place between A and C). The time unit is . Note the duration of the preparation process with method 2, and the difference in final equilibrated states compared to method 1.

The second method to prepare numerical samples allows for aggregate formation before imposing an external pressure. Along with method 1, its different stages are schematically presented on Fig. 4.

The aggregation phenomenon plays an important role in the experimental preparation procedure of refs. Valverde et al. (2001a, b) in which powder particles in a fluidized bed collide and stick to each other. Then they settle under their weight when the upwards air flow is abruptly shut off. The numerical method was designed to reproduce, in some idealized way, the final state of a set of colliding particles in the absence of external force fields. In the initial disordered low-density configuration (the same -state as in method 1), particles are now attributed random velocities drawn according to a Maxwell distribution, with mean quadratic velocity .

We performed systematic sets of simulations of disk packings with (see Eqn. (7)). is thus large enough for the initial kinetic energy to overcome potential energy barriers in the process of aggregation. (The dependence of the final packing structure on this initial velocity of agitation, or “granular temperature”, in the assembling stage in systems with small RR will be studied in Section IV.2.6).

Once launched with such random velocities the particles are left to interact and stick to one another within a cell of constant size, forming larger and larger aggregates, as appears on the image marked “B2” on Fig. 4. Eventually, all particles are connected to one another by adhesive contacts, and reach an equilibrium position. At this stage, the two degrees of freedom of the cell are set free, and the stress-controlled calculation proceeds with (or ) until an equilibrium state is reached. This relaxation step does not lead to any rearrangement of the contact structure, it only entails a very small increase of the solid fraction (hence the values slightly larger than given below). The final equilibrium structure exhibits large density inhomogeneities, as apparent on Fig. 4, which are characteristic of aggregation processes Smirnov (1990), and will be quantitatively studied in Section V.

Unlike cohesionless systems, which are devoid of any “natural” state of stress, clusters of cohesive particles can exist in a well defined state of mechanical equilibrium in the absence of any external force. Once the state at zero pressure is obtained, we subsequently apply the same load as in method 1, which results in further compression and notable changes in the packing structure: increases from values close to up to the range (see fig. 4). Nevertheless, the final solid fraction under is considerably lower than the one obtained with method 1.

It should be noted on Fig. 4, which summarizes the assembling procedures, that the aggregation stage makes method 2 computationally quite costly because of the time necessary for clusters to merge, and especially for the stabilization of loose samples in equilibrium configurations (lower contact numbers implying lower rates of energy loss as well as larger and slower fluctuations of soft, tenuous structures). In an attempt to limit the influence of compaction dynamics, which results in denser samples when the lower density of the initial state allows the compaction process to accelerate more (as noted in Sec. III.1), we tested the effect of limiting the maximum strain rate . Without any limitation, we obtained a maximum value . Using the samples with (the lowest value used in this work) with , three different values for were tested: , and . The condition gave a final state close to the original one. The others two values produced similar results, with a relative decrease in density of about compared to the original procedure. We chose to enforce condition , to save computational time. This value has been applied to prepare all samples studied in the following.

Fig. 4 shows that method 2 succeeds in stabilizing open structures. Final solid fractions agree with the typical values observed in powders if one uses the correspondence between 2D and 3D packing fractions suggested by Campbell in Campbell and Brennen (1985):


Numerical samples under , with solid fractions around , would correspond to a powder consolidated in the laboratory under 1 Pa with a solid fraction of about . This is in satisfactory agreement with the experimental results of Ref. Valverde et al. (2004).

We therefore regard method 2 as an appropriate way to reach an essential objective of this work, since stable loose structures are obtained.

Although we perform simulations of a mechanical model, the final configurations exhibit at first sight (Fig. 4) similar features as those obtained with geometric algorithms implemented in numerical studies of colloid aggregation models Meakin (1999); Kadau et al. (2003). We are not aware of similar results in the literature, at least with equilibrium requirements comparable to those of Sec. II.4.

Tenuous, fractal-like contact networks contain denser regions and large cavities. Such heterogeneities produce long range density correlations, to be analyzed in Sec. V. Without tensile contact forces, the walls of the cavities, comprising particles that are pushed towards the hole by the resultant of contact normal forces, would tend to buckle in.

We regard method 2 as yielding typical results for assembling processes in which particles form tenuous aggregates before they are packed in a structure that is able to sustain a confining stress. In the sequel, we focus on the tenuous structures obtained with method 2.

Global characterization of loose packings at and .

We simulated four samples with 1400 disks and three of 5600 for , rather than lower initial densities, in order to achieve statistical significance at affordable computational costs, and to check for possible size effects. This set of samples will be denoted as series A.

Samples with (series A0), which require the initial cell to shrink more before a stable network can resist the pressure, request longer calculations. Although some samples were prepared at , we do not use them any more in the following, except for the values showed in Table 4.

To accelerate the numerical assembling procedure, we also created samples with , using an intermediate value of , and softer contacts, such that in the initial aggregation stage (recall the time step is proportional to ). Once equilibrium was reached with , we slowly changed the stiffness parameter from to , and recorded the final equilibrated configuration. This procedure is about ten times as fast as the normal one, and generates similar structures and coordination numbers as series A prepared with the same -state density. We shall refer to this set as series B.

In table 4 we list the corresponding results for solid fractions and coordination numbers. In such data we did not find a significant difference between the two different sample sizes, and therefore we did not distinguish between sizes in the presentation of statistical results.

no RR RR
(series A0)
(series B)
(series A)
Table 4: Values of and on equilibrating configurations at with . Samples with correspond to , the others to .

The tenuous networks obtained with method 2 collapse on changing the pressure: Table 5 gives the new values of and after the compaction caused by the pressure increase from to .

no RR RR
(series B)
(series A)
Table 5: Values of and in equilibrated configurations at . These results are averaged over the whole set of samples prepared with , for (series B) on the one hand, and for (series A) on the other hand. Series A0, prepared with , yielded very similar results but due to computational costs the number of samples was too small to record data in statistical form.

Structural changes between and are shown on Fig. 5, which illustrates by means of four selected snapshots the mechanism of the closing of pores in a 1400 disks sample of series A. The first image corresponds to equilibrium at , and the fourth one to equilibrium under . The two others show intermediate, out of equilibrium configurations during the collapse. One may appreciate how the denser regions grow and merge while pores shrink.

Figure 5: Configuration of 1400 disk sample of series A without rolling resistance. Note the gradual closing of pores as the external pressure is increased from (first image) to (last image) going through two intermediate stages.

Fig. 5 also makes it quite evident that the size of 1400 disk samples is not very much larger than the scale of density heterogeneities (typical diameter of large pores or dense regions, which will be studied in Sec. V). These systems will exhibit large fluctuations in their mechanical properties: the rectangular shape of the final configuration displayed on Fig. 5 shows that the disorder is large enough for the mechanical response of the system to become anisotropic. Isotropy should be recovered in the limit of large sample sizes, .

Finally, Fig. 6 displays the histogram of local coordination numbers (percentage of particles interacting with others, ), for the same samples as those of Tables 4 and 5 (, ).

Figure 6: Distribution of local coordination numbers in loose samples of series A obtained with method 2. Samples of series B gave a similar distribution.

It is remarkable that this distribution, in spite of the large difference in sample geometries, remains rather close to the one observed in the denser packings made with method 1 (compare , case), just like global coordination numbers take very similar values in samples prepared with both methods (see Tables 3 and 5), in spite of the very different solid fractions.

An essential conclusion of the present study is therefore, for one given material, the absence of a general relation between the density of a cohesive packing and its coordination number, in spite of previous claims Yang et al. (2003a). Both quantities are determined, rather, by the conjunction of micromechanical laws and sample preparation history.

Effects of micromechanical parameters

Adhesion should enhance the role of sliding friction and rolling friction, because the limiting values for tangential contact forces and rolling moments are both proportional to the elastic repulsive part of the normal force, ( , ). Consequently, contacts with the equilibrium value of the elastic deflection for an isolated pair of grains transmit no normal force, but are able to sustain tangential force components as large as and rolling moments as large as . Those values might turn out to be large in comparison to the typical level of intergranular forces under low external pressure ().

Figure 7: Typical configurations of 1400 disk samples of series A with (left) and without (right) rolling resistance, at (a) and (b). Note the difference in local structure of thin “beams” joining dense regions with or without RR.

Therefore, even very low values of and should affect the final structure of equilibrated packings considerably more than in the cohesionless case. This is indeed the case for the coordination numbers observed in our simulations (see tables 2 and 3) which dropped more significantly, upon introducing the small level of RR we have been using, in cohesive systems than in cohesionless ones.

On Fig. 7 we show the configurations at (a) and (b) of the same sample assembled using method 2 with RR (left) and without RR (right). The denser regions in the inhomogeneous packings are joined by slender “arms” (see Figs. 5 and 7). Such arms can in principle reduce to a chain of particles in the presence of rolling resistance. Such chains are otherwise destabilized by a rolling mechanism, hence the difference in the thickness of the arms with or without RR (see the blown-up detail in Fig. 7-a), the lower coordination numbers of configurations assembled with RR (Tables 4 and 5). This might also explain the greater fragility of equilibrium configurations with RR, in which a larger compaction step (see Table 5) is necessary, on applying , before a new stable structure is reached.

Figure 8: Final coordination number versus initial quadratic average velocity in agitation stage of method 2, normalized by characteristic velocity . The arrow points to the value most often used in our calculations.

Another important parameter is the initial velocity of agitation, . Its influence has been assessed on one 1400 disks sample, with . The changes of coordination number with at are presented on Fig. 8.

Low velocity values produce more tenuous aggregates (), since even a small level of RR is able to slow down local rearrangements and stabilize tree-like structures (i.e., devoid of flops) immediately after the collisions between particles or small clusters.

A large kinetic energy cannot be absorbed by the RR, and as a result disks are able to rotate, which leads to better connected structures (). In a sense, a large kills the effects of RR, and packings are similar to those made without RR in such cases.

We therefore conclude that the connectivity of loose samples with RR assembled by aggregation depends on the initial magnitude of velocity fluctuations and on the level of rolling friction.

As figure 8 shows, the same trend was found on reducing contact stiffness parameter , as a larger translational and rotational compliance creates more contacts.

is analogous to the particle fluctuating velocity in experiments on gas-fluidized beds of xerographic toners under gravity Valverde et al. (2001a). Such velocities are larger than the gas velocity by two orders of magnitude. Typically, one has mm/s, while , deduced from the contact parameters with relation (7) is about 1 cm/s. Such a value is therefore comparable to the particle fluctuation velocity.

Of course, such a comparison is only indicative, because the influence of on packing structures depends on , and is also very likely to be affected to some extent by the viscous dissipation model we have adopted. Both rolling resistance and viscous forces are micromechanical features for which no accurate physical identification is available. Yet, it seems plausible that powder packings, because of their initial agitated states, stabilize in better connected states than predicted by geometric aggregation models.

We now turn our attention in the next section to the forces in the contact networks, in particular the loose ones formed with method 2.

Iv Mechanical characterization of contact networks

Many numerical studies, in the past 15 years, have addressed the issue of contact network geometry and force distribution in cohesionless systems Radjai et al. (1996). The image of force chains, i.e. a pattern in which larger intergranular forces tend to line up on the scale of several grains, was evidenced in experiments Mueth et al. (1998); Blair et al. (2001) and simulations Ouaguenouni and Roux (1997); Makse et al. (2000), and the p.d.f of contact force values has often been measured and studied. An interpretation of the mechanical role of “force chains” Radjai et al. (1998) is that they carry the essential part of deviatoric stress, while the contacts carrying the lower forces are less sensitive to stress orientation and laterally stabilize the strong force chains against buckling.

The main features of the distribution of forces and their spatial correlations have been reproduced by approximate models Coppersmith et al. (1996) based on local equilibrium rules on each grain, supplemented by inequality constraints. One important such constraint is released in cohesive systems, in which normal force components can have either sign. It is therefore worth investigating how the usual features of force-carrying structures in equilibrated granular packings are affected by the presence of negative normal forces. One may also wonder to what extent the considerable difference in the density fields will affect the force patterns, given that the coordination of the force networks, as observed previously, does not seem to be very sensitive to density levels and density fluctuations.

iv.1 Force scale and force distribution

The first obvious distinction between cohesive and non-cohesive systems is the appearance of a new force scale , in addition to the one provided by the confining pressure, i.e., , the ratio of those characteristic forces defining the reduced pressure, . It is especially interesting to investigate the values and spatial organization of forces in systems with , as little information is to be found in the literature on this issue: numerical studies of loose cohesive systems Wolf et al. (2005) tend to focus on density and geometry of packings as a function of applied stresses. Some information on force networks is provided in a recent publication Richefeu et al. (2006) on bead assemblies with capillary cohesion, but the confining stress is considerably higher is that study ( of order 1) than in the present one.

In the absence of cohesion, the distribution of force values is usually presented in a form normalized by its average, which itself scales with the applied pressure. This scaling can be made more quantitative on using a general relation between pressure and the average normal contact force and particle diameter , which is known in the literature on powders as the Rumpf formula. We write it here in a form involving the spatial dimension , which is valid both for and :


In (15), stands for the typical grain diameter. This relation can be made more accurate if one notes that it stems from the standard formula for stresses in an equilibrium configuration (see the r.h.s. term in Eqn. 3). To derive the formula, defining , the average pressure, one assumes and then neglects correlations between particle radii and forces, assuming


Then, with a simple transformation of the sum, one obtains


With and our diameter distribution (for which and ) this yields


We found relation (18) to be remarkably accurate in all our simulations, with or without cohesion, with configurations obtained by either method 1 or method 2, thereby checking that the correlations between particle sizes and contact forces could safely be neglected on writing (16).

Without cohesion, Eqn. (18) yields the correct scale for forces, i.e. the frequency of occurrence of intergranular forces larger than a few times is very small. With cohesion, when or , contact forces of order are quite common, as shown on Fig. c, on which normal force distributions are represented.

(a) No cohesion,
(b) With cohesion, method 1, , and
(c) With cohesion, method 2, and ,
Figure 9: Distribution of normal forces for series A samples. The non-cohesive case (a) is normalized by the average repulsive elastic part. The cohesive cases (b and c) are normalized by (note that the average of the elastic part of is in cohesive cases with ).

Hence Eqn. (18) cannot be used to predict “typical” contact forces. The presence of forces of order explains the sensitivity of type 1 and type 2 samples with to friction coefficient and rolling resistance: densities and coordination numbers (tables 3, 4 and 5), in cohesive systems prepared under or with and with , or with and without RR, differ significantly. Otherwise, if contact forces were of order of the average , the value of which is correctly predicted by (18), thresholds or even would be very large compared to typical forces and moments, and become irrelevant.

(It should be recalled that Rumpf’s name is often associated (as in Ref. Richefeu et al. (2006)) to a means to predict the macroscopic tensile strength of a powder. As the essential ingredient of the Rumpf approach Rumpf (1958) is Eqn. 15, we refer here to that relation (like in Castellanos (2005)), as the Rumpf formula).

Normal force distributions in cohesionless, cohesive type 1 and cohesive type 2 samples, the latter being obtained with (series A), are shown on Fig. c. Those distribution functions are roughly symmetric about , decay approximately exponentially at intermediate values, and vanish at , and . In type 2 samples without RR, for , there is a finite proportion of contacts carrying vanishing forces, about one fourth in the A series (). In addition to this Dirac mass, there might be a power-law divergence near 0, with an exponent our level of statistics is not sufficient to resolve accurately (about 0.6 to 0.8 in the range of forces between and ). This proportion of zero forces is smaller, down to , with RR, and drops as reaches , to and , respectively, without and with RR. It is worth pointing out that the corresponding contacts carry zero total forces, i.e. vanishing normal components (, see (11)) and no tangential elastic displacement either. In principle we cannot distinguish them from forces below the numerical tolerance defined in Sec. II.4. However, as we shall argue below in Section IV.2, under one could expect all contact forces to vanish, and non-zero forces are related to the small, but finite degree of force indeterminacy.

Before turning our attention to such features and to the spatial organization of forces, let us briefly discuss the differences between sets of (type 2) samples A and B. B samples, which are obtained with the “accelerated” procedure and (see Sec. III), exhibit, due to their specific history, larger forces at , with as many as of the contacts transmitting normal forces such that , while this proportion lies below in A samples. On the other hand, B samples are looser, with more open contact networks under and a larger proportion of contacts (about one third in configurations without RR) carrying vanishing forces. In the following we shall use them to illustrate qualitative tendencies in very loose samples.

When the pressure is increased to , differences in force distributions between A and B samples, despite their different solid fractions (see table 5), have considerably decreased, as shown on Fig. LABEL:fig:hist2636.

The influence of such differences in the aggregation stage as those between our samples A and B are therefore expected to fade out after the systems are compressed to higher pressures and densities.

iv.2 Packing structure and force patterns

The spatial organization of forces in type 2 samples, which we now discuss, is related to the distribution of force values, and should determine the ability of given configurations and contact networks to support stress increments. We first discuss systems without RR, then with the small RR values we adopted in most cases (see Table 1). We emphasize the role of force indeterminacy and assembling history (the collisions by which cohesive clusters were built) in the final force patterns in equilibrium under vanishing or low applied stress. Extreme cases of systems with large RR on the one hand, or without friction on the other hand, are useful reference situations, which we briefly examine and discuss. We conclude this part with a discussion of the main physical implications of the relationships between force patterns, assembling process, geometry and micromechanical parameters

Qualitative aspects of force networks with no RR

It is instructive to represent the forces carried by the contact network with a visualization of positive (repulsive) and negative (attractive) normal forces, as was done on Fig. b, showing the force network in one type 1 sample. Figs 11 and 12

Figure 11: (Color online) Sample of type 2 (N=1400), in equilibrium under after aggregation stage, with solid fraction (series B). Same conventions as on Fig. b, except for the blue color corresponding to contacts carrying a total force below tolerance (deflection and no mobilization of tangential force). Note the large number of such interactions and the local compensation of attractions and repulsions in small prestressed clusters. To help visualize unstressed regions, disks only interacting at contacts bearing forces below tolerance are filled in light grey.

respectively correspond to equilibrated samples prepared with method 2 under (immediately after the aggregation stage) and under , without RR. They are represented here with (approximately) the same scale. Both belong to the () B series. Line widths, which are proportional to the intensity of the total interaction force, i.e. to , witness the presence, in spite of the low pressure, of many forces of order (which correspond on the figures to line thicknesses comparable to particle radii). Stressed clusters, in loose type 2 samples under , are separated by large parts of the interacting network in which contacts carry vanishing forces: the corresponding normal deflection is (Eqn. (11)) and there has been no elastic relative tangential displacement.

Figure 12: (Color online) Sample of Fig. 11, with same scale and color conventions, in equilibrium under . The solid fraction increased to . The threshold force (used to distinguish blue lines and grey disks) was set to .

Attractions (green) and compressions (red) have to compensate for Eqn. (18) to hold true. This compensation appears to operate on a smaller scale in type 2 samples, because internal forces were previously balanced within isolated particle clusters. Such a local balance of forces is quite conspicuous at (Fig. 11), in which internal stresses in small clusters often take the form of a peripheral tension compensating a radial compression, or the other way round. This contrasts with samples prepared with method 1 (Fig. b), in which the spatial distribution of forces is more similar to the familiar “force chain” pattern of cohesionless systems, although there are of course compressive and tensile “chains”. Unstressed regions are rather scarce in type 1 samples, although some areas with smaller forces are still present. The structure of type 2 samples under (Fig. 12) is somewhat intermediate: isolated stressed clusters are still present, but elongated, force-chain-like structures emerge.

To characterize such force patterns in a slightly more quantitative way, one can evaluate a threshold force such that contacts carrying a force with percolate through the sample. Such a criterion was used to identify a “strong” subnetwork of force chains in Radjai et al. (1998). One observes of the order of the tolerance in samples with no RR and , which shows that stressed regions are isolated “islands” within the network. raises to slightly less than under .

Configurations of series A, assembled with , possess the same qualitative features, although quantitatively slightly weaker, due to their higher density. For instance, local stressed regions are somewhat less isolated, with a threshold force between and at .

Force indeterminacy (without RR)

The presence of large interaction forces of order in equilibrated samples is not obviously necessary, and is related to the assembling process. Let us imagine particles are brought very slowly, one by one, within interaction range of the previous network, thus gradually building a unique cluster in equilibrium in the absence of external stress. One could expect, rather, each new contact to stabilize with and . The existence of non-zero interaction forces in equilibrium is related to the hyperstaticity or force indeterminacy of the contact network. On writing all equilibrium equations for grains and collective degrees of freedom (i.e., setting acceleration terms to zero in Eqns. 1, 2 and 3) and regarding all contact forces as unknowns, the degree of force indeterminacy h (or degree of hyperstaticity) is the number of remaining independent unknowns, which cannot be determined by the equilibrium requirement. If , knowing that some equilibrium forces exist (since an equilibrium state has been found), then one would necessarily have all interaction forces equal to zero under (since this is one obvious possible solution). The notion of force indeterminacy has been recently discussed by different groups in the context of granular materials, essentially because of the special case of rigid frictionless grains, for which the contact network is generically such that forces are uniquely determined Ouaguenouni and Roux (1997); Roux (1997); Moukarzel (1998); Tkachenko and Witten (1999); Roux (2000). The degree of force indeterminacy is linked to the number of degrees of freedom, equal to (or if the cell sizes can change), to the number of contacts , the number of distant interactions and the number of independent mechanisms or floppy modes k (also called degree of hypostaticity Roux (2000)) by the following relation (written here for a fixed cell)


A proof of this simple result (which is classical in structural engineering), and the relation of numbers h and k to the rigidity and stiffness matrices of the contact network, are recalled in Appendix A. Mechanisms are those sets of velocities (or small displacements, dealt with as infinitesimal) which entail no relative velocities (or small relative displacements) in contacts. For distant interactions, only normal relative velocities are relevant, hence their particular treatment in (19). In Appendix A we explain how we determine whether a given configuration is rigid, i.e., devoid of mechanisms (apart from the two global translational motions of the whole set of grains, rendered possible by the periodic boundary conditions). It is customary to relate the level h of force indeterminacy to the coordination number in granular materials. However, this is not possible in general, which motivated our recalling (19) in its complete form. (19) can be rewritten, neglecting the very scarce distant interactions, as

Hence, in the absence of floppy modes, . However, there are still a few floppy modes on structures like those of Fig. 11 at the end of the aggregation stage, and this relation, which predicts a small degree of hyperstaticity relative to the number of grains (see Tables 4 and 5), is only approximate. Some mechanisms are due to the (exceptional) 1-coordinated disks and others, less trivial, are associated with larger parts of the structure which are connected to the rest of the packing via one single 2-coordinated disk. This floppiness is obviously related to the assembling process: before any external pressure is applied, nothing really requests the aggregates to possess a rigid backbone. The free motion of mechanisms in assembling method 2 is largely responsible for the very long equilibration time (see Fig. 4): such motions entail no restoring force and no dissipation of kinetic energy. Floppy modes in the final state obtained with our criteria (Section II.4) being scarce (typically, a few such mechanisms per 1400-disk sample), we conjecture that they would disappear entirely on adopting stricter equilibrium requirements in terms of kinetic energy. If a mechanism survives, it should generically be in motion with a non-vanishing velocity, as a residual effect of the initially agitated state. As the connected aggregate partly folds onto itself, such motions should eventually create new adhesive contacts, thereby reducing k, until the network becomes rigid. Once some rigid aggregate is formed in the assembling process, it will keep the same shape and structure, unless the collisions and perturbations it subsequently undergoes cause it to break, because of the limited tensile strength of contacts or because of the Coulomb inequality. This is the reason why the initial mean quadratic velocity of isolated grains in method 2 should be compared to , as given by (7).

It is easy to see that the closing of one contact can convert an aggregate from floppy to hyperstatic, the simplest example thereof being the “double triangle” structure of Fig. 13. By (19), this small structure, which is rigid ( counting the free motions of an isolated object in 2D) has a degree of force indeterminacy .

Figure 13: Three (grey) disks initially forming an isostatic structure, when a fourth one (coming from the left) adheres to one of them (dotted position), it can roll (this is a floppy mode) until another (fifth) contact is formed, stabilizing it in the final position drawn with continuous line. The final “double triangle” structure is hyperstatic. Final forces (see main text) were computed for different initial velocities of the mobile disk and for different values of impact parameter .

This is how the self-stressed clusters of Fig. 11 are formed. Such structures have a strong influence on force values and force distribution. In particular, we show now that they entail specific correlations between normal and tangential force components in contacts.

Local patterns and specific force orientations

Fig. 14 shows all contact force values as points in the plane for a 5600 disk sample without RR equilibrated under .

Figure 14: Values of normal and tangential contact forces in a 5600 disk, type 2 sample, in equilibrium under , with (A configuration). In addition to the remarkable cross-shaped pattern, marked with dashed lines of slopes , note the large number of very small forces, the numerous points with and the relevance of the value of the friction coefficient ( here), as a small number of forces approach the Coulomb cone.

Fig. 14 displays a striking X-shaped distribution in the plane, corresponding to a ratio of This cross pattern fades away in systems which have rearranged to support , although the corresponding specific ratios are still overrepresented, as shown on Fig. b.

Figure 15: Histograms of angle , between normal vector and total contact force . Conventionally, for a repulsive normal force and , and for a tensile normal force and . Shaded histograms (grey) correspond to B configurations (), bold-line non-shaded ones (black) to A ones ()

The cross pattern of Fig. 14 corresponds to angles and on Fig. b, and the second graph shows that still corresponds to a peak in the distribution once the sample has been compressed (and rearranged) to . As other characteristic features of force patterns in loose type 2 samples, this correlation between tangential and normal force components is stronger in the more tenuous networks of series B samples, for which the data are also represented on Fig. b. The difference between both sample series tends to disappear on compressing to (Fig. b).

The prevalence of ratio is in fact easy to understand. Many disks are in equilibrium with two contact forces, with two other disks which are themselves contacting each other, as on Fig. 16.

Figure 16: The bottom disk, marked D, of radius , is in contact with two other disks 1 and 2, themselves touching, whose radii are and . At equilibrium, contact forces on disk D should be carried by the dotted line joining its 2 contact points, which determines the ratio of tangential to normal force components.

(Occasionally, a third contact might be present, bearing a much smaller force, which we neglect in the present argument). In such a situation, without RR, the three equations expressing the balance of forces and moments on disk D involve four unknown force components. Labels corresponding to contacts with disks marked 1 and 2 like on Fig. 16, one obtains, on counting positively repulsive normal forces on disk D and tangential forces with a positive moment:


The ratio in (20) varies for the radius distribution we are using in the present study but its most frequent value, corresponding to , is . Fig. 17 shows the same graph as that of Fig. 14, in the case of a loose packing of disks with the same radii. In agreement with formula (20), the “X” shape is sharply defined.

Figure 17: Values of normal and tangential contact forces in a 5600 disk, type 2 sample of monodisperse disks in equilibrium under . Note the sharp“X” shape on blown-up detail of small forces.

To understand the frequency of occurrence of very small values, let us now consider again the smallest cluster with force indeterminacy, without RR, which comprises four disks and 5 contacts, as schematized on Fig. 18. Fig. 18 shows graphically that the balance of contact forces implies that the tangential force within the contact corresponding to the common base of the two triangles should be very small, thereby explaining the “dense line” along the axis on Fig. 14. It can be checked by direct inspection that local simple patterns like those of Figs. 16 and 18 are indeed typical for the forces with ratios around , or with .

Figure 18: Hyperstatic 4-disk cluster, with 5 contacts. The force at the contact point between 1 and 2 should be carried by the continuous line joining this point to the intersection of the dotted lines. Those lines are respectively defined, as on Fig. 16, by the two contact points