Computer simulation of model cohesive powders: Plastic consolidation, structural changes and elasticity under isotropic loads
aa r X i v : . [ c ond - m a t . m t r l - s c i ] O c t Computer simulation of model cohesive powders:Plastic consolidation, structural changes and elasticity under isotropic loads.
F. A. Gilabert, ∗ J.-N. Roux, and A. Castellanos Faculty of Physics, University of Seville, Avda. Reina Mercedes s/n, 41012 Seville, Spain. Universit´e Paris-Est, Institut Navier,Laboratoire des Mat´eriaux et des Structures du G´enie Civil † ,2 All´ee Kepler, Cit´e Descartes, 77420 Champs-sur-Marne, France. (Dated: October 27, 2018)The quasistatic behavior of a simple 2D model of a cohesive powder under isotropic loads isinvestigated by Discrete Element simulations. We ignore contact plasticity and focus on the effectof geometry and collective rearrangements on the material behavior. The loose packing states, asassembled and characterized in a previous numerical study [Gilabert, Roux and Castellanos, Phys.Rev. E , 011303 (2007)], are observed, under growing confining pressure P , to undergo importantstructural changes, while solid fraction Φ irreversibly increases (typically, from 0.4–0.5 to 0.75–0.8). The system state goes through three stages, with different forms of the plastic consolidationcurve , i.e. , Φ as a function of the growing reduced pressure P ∗ = P a/F , defined with adhesionforce F and grain diameter a . In the low-confinement regime (I), the system undergoes negligibleplastic compaction, and its structure is influenced by the assembling process. In regime II thematerial state is independent of initial conditions, and the void ratio varies linearly with log P [i. e.∆(1 / Φ) = λ ∆(log P ∗ )], as described in the engineering literature. Plasticity index λ is reduced inthe presence of a small rolling resistance (RR). In the last stage of compaction (III), Φ approachesan asymptotic, maximum solid fraction Φ max , as a power law, Φ max − Φ ∝ ( P ∗ ) − α , with α ≃ ξ of fractal density correlations decreases, force patterns reorganize from self-balanced clustersto force chains, with correlative evolutions of force distributions, and elastic moduli increase by alarge amount. Plastic deformation events correspond to very small changes in the network topology,while the denser regions tend to move like rigid bodies. Elastic properties are dominated by thebending of thin junctions in loose systems. For growing RR those tend to reduce to particle chains,the folding of which, rather than tensile ruptures, controls plastic compaction. PACS numbers: 45.70.-n,81.40.Lm,61.43.Hv,83.10.Rs
I. INTRODUCTION
Cohesive granular materials are present in many nat-ural or industrial processes, the understanding of whichrequires studies of their rheology under small confiningpressures, when tensile intergranular forces play a majorrole. In such cases cohesive materials exhibit specific fea-tures that do not exist in cohesionless grain assemblies,such as the ability to form stable structures at low den-sity and the sensitivity to stress intensity , as opposed tostress direction. Macroscopic constitutive laws and phe-nomenological tools have been developed and used in sev-eral engineering fields: mechanics of cohesive soils (claysand silts) [1, 2, 3, 4], metallic powder processing [5],modeling and treatment of ceramic powders [6, 7, 8, 9],handling of xerographic toners [10]. One simple mate-rial is the assembly of wet beads [11, 12, 13], in whichsome microscopic observations are possible [12, 13]. How-ever, wet grain packs are only slightly less dense than † LMSGC is a joint laboratory depending on Laboratoire Centraldes Ponts et Chauss´ees, ´Ecole Nationale des Ponts et Chauss´eesand Centre National de la Recherche Scientifique. ∗ Electronic address: [email protected] dry ones, and do not enable the study of loose structuresobtained with powders. In general, the behavior of ma-terials under proportional load (oedometric or isotropiccompression) is characterized by the consolidation curve,which describes the irreversible compaction under grow-ing stress [1]. Density can increase by factors of 3 or 4under growing load.Although numerical simulations have been widely usedfor several decades [14] to investigate microscopic mech-anisms and classify mechanical properties of granularsystems, studies of cohesive materials are still far lesscommon, and almost exclusively limited to dense ma-terials. Thus, the effects of capillary cohesion in wetsand or bead packs have been simulated [15, 16], aswell as the compaction of ceramic and metallic pow-ders [17, 18, 19, 20, 21, 22, 23] to states of very highdensity, or the behavior in shear tests of 2D dense cohe-sive packs with plastic deformation of contacts [24, 25].Loose structures formed by particles packed under grav-ity and stabilized thanks to adhesion have been simu-lated [26]. Of particular relevance to the present study,among the very scarce numerical studies of loose pack-ings [27] stabilized by cohesion and of their collapsingunder growing loads, are the works by Bartels, Kadau,Wolf et al. [28, 29, 30, 31] on the oedometric compres-sion of granular assemblies with initial low densities. Thisresearch group studied a dynamical compression regime,and observed a shock wave propagating through the sam-ple. Shear flows of cohesive granular materials have alsobeen simulated [32, 33, 34, 35].In a previous article [36], hereafter referred to as pa-per I, we studied by numerical simulation the assemblingprocess, the structure and the force patterns of a model,two-dimensional (2D) cohesive granular material in looseequilibrium configurations. We now investigate the me-chanical behavior of the same model granular materialin isotropic compression and pressure cycles, as well asthe evolution of various characteristics of intermediateequilibrium states as plastic compaction proceeds.As in paper I, we keep the external pressure as themain control parameter. The adhesive strength F incontacts sets a force scale in the material behavior, andhence (in 2D) the reduced pressure, defined as P ∗ = aPF , (1)in which a is a typical grain diameter, is a crucial di-mensionless state parameter. The main objective of thepresent paper is the study of the process by which, aspressure is increased, cohesion-dominated loose struc-tures, for which P ∗ ≪
1, get irreversibly compactedas P ∗ increases until pressure dominates ( P ∗ ≫ e.g., inRef. [31]. However, our approach produces homogeneous,isotropic, equilibrium configurations under varying loadand is therefore apt to provide more detailed informationabout the connections between macroscopic constitutivelaws and microstructural or micromechanical features.The present paper is self-contained and can be under-stood without reading paper I. A summarized descriptionof the material properties and of the initial configura-tions (studied in paper I) is provided in Section II. Themacroscopic material response in isotropic compression,with the possible influence of the initial state proper-ties, is studied in Section III. Then, various microscopicaspects of the consolidation process are investigated inthe sequel: density correlations (with their fractal behav-ior over some length scale [36]) are investigated in Sec-tion IV, force networks and force distributions are dealtwith in Section V, while Section VI focuses on elasticmoduli. Section VII discusses qualitatively some micro-scopic aspects of the consolidation behavior. The finalsection, part VIII, summarizes the results and suggestsdirections for future work. Sections IV and V can, atfirst, be read independently from each other. The sameremark applies to Sections VI and VII. II. MODEL MATERIAL AND SIMULATIONPROCEDURESA. Definitions and basic equations
The material and the simulation method are identicalto those of paper I [36], which the reader might refer to foradditional technical details, and for a physical discussionof some of the model ingredients. For the sake of com-pleteness, we however provide a summarized descriptionbelow. The contact law is an elaboration of the oftenemployed spring-dashpot model with Coulomb friction,in which two additional ingredients are introduced: anattractive force and, possibly, some resistance to rollingat contacts. The model material is a 2D assembly ofdisks with diameters uniformly distributed between a/ a , enclosed in a rectangular cell with periodic bound-ary conditions in both directions. Both lengths L , L defining the cell size and shape are variable, and sat-isfy equations of motion designed to impose given valuesof diagonal stress components σ = σ = P . Stressesare controlled by a variant of the Parrinello-Rahmanmethod [37]. In equilibrium, both diagonal stress com-ponents σ α , ( α = 1 , A is the sample surface area): σ α = 1 A X ≤ i Grains interact with forces of elastic, adhesive, fric-tional and viscous origins. The static part of the normalcomponent F ijN of the force transmitted by grain i to itsneighbor j is a function of h ij , the distance separatingdisk perimeters. A negative h ij means that the grainsoverlap, in which case they repel each other with a nor-mal elastic force F e,ijN = − K N h ij . This force vanisheswhenever h ij > 0. (Overlap h ij < F a,ijN , equal to − F for contacting disks( h ij < F a,ijN has a finite range D , fixed to 10 − a ,and varies linearly between − F and zero as h ij growsfrom 0 to D . F is the maximum tensile force a contactmight support without breaking off. The normal contactlaw thus introduces a force scale, and a dimensionless pa-rameter, the stiffness parameter , κ ≡ aK N /F . κ charac-terizes the amount of elastic deflection h under contactforce F , relative to grain size a ( h /a = κ − ). κ is setto a large value, κ = 10 , so that the elastic deflectionsin contacts remain so small that they can be neglected incomparison to all other length scales in the problem (in-cluding interstices between neighbors [38]). The packinggeometry can be regarded as that of an assembly of rigidgrains (as formally dealt with in the “contact dynamics”simulation method used in [31]).To the static contributions F eN and F aN to the normalforce we add a viscous damping term opposing the rel-ative normal velocity of i and j when the disks touch( h ij < e N in binary collisions if F isset to zero. e N is set to a low value, e N = 0 . 015 inour simulations. In the presence of attractive forces theapparent restitution coefficient in a collision will dependon the initial relative velocity. For small kinetic ener-gies the particles will eventually stick to each other. Theminimum receding velocity for two particles of unit mass(the unit mass is chosen equal to the mass of a disk ofdiameter a ) to separate is V ∗ √ 2, with V ∗ = p F D . (3)The elastic tangential force in contact i, j , F ijT , is to beevaluated incrementally. In case of no tangential sliding,it varies linearly with the relative tangential displacementat the contact point, involving a tangential stiffness con-stant, K T . In the case of sliding, which occurs when theelastic law would cause F ijT to pass one of the Coulombbounds ± µF e,ijN , then F ijT stays equal to ± F e,ijN . Therelative tangential displacement at the contact point in-volves displacements of disk centers and rotations. TheCoulomb condition introduces the friction coefficient , µ .It should be pointed out that it applies to the elasticrepulsive part of the normal force only. Thus, a pairof contacting grains with h ij equal to F /K N = h , theequilibrium distance, such that the sum of elastic and ad-hesive terms vanishes, can transmit a tangential force F T such that | F T | ≤ µF . (The importance of this feature ofthe contact law for collective properties macroscopic be-havior of particle assemblies was stressed in paper I forisotropic, static states, and in Ref. [34] in steady-stateshear flows). All simulations reported here were carriedout with µ = 0 . rolling resistance (RR) atcontacts, which is modeled as in [39]. Two additionalparameters are necessary: a rolling spring constant , K R ,with dimension of a moment, expressing proportional-ity between relative rotation and rolling moment (i. e.,a torque concentrated at the contact point), as long as the rolling friction threshold is not reached; and a rollingfriction coefficient , µ R with the dimension of a length,setting the maximum absolute value of the rolling mo-ment Γ R to µ R F eN , proportional to the elastic part of thenormal force. The implementation of this rolling law isanalogous to that of the tangential one, with the rollingmoment and the relative rotation respectively replacingthe tangential force and the relative tangential displace-ment. A contact for which the total normal force is equalto zero in equilibrium, with F eN = K N h = F , maytransmit a rolling moment Γ R with | Γ R | ≤ µ R F eN . Sincepoint contacts do not transmit torques, the rolling resis-tance stems from the irregularity of grain surface. Twocontacting grains touch each other, in general, by twopoints (in 2D), which are separated by some microscopicdistance l that is characteristic of the particle shape. µ R should be proportional to l , and K R proportional to l .We set µ R = µl and K R = K N l , with, in most calcula-tions with RR, l = a/ µ e N κ K T K N D a K R K N a µ R a . . 015 10 − − . tions were also performed with larger RR (up to l = a , µ R = 0 . a ). C. Initial states In paper I, two extreme cases were studied in the as-sembling stage of cohesive packings under low P ∗ . First,an N -particle sample of hard-disk fluid is prepared atsolid fraction Φ I in a fixed cell. Then, in type 1 sys-tems, velocities are set to zero and the external pres-sure control is started, until an equilibrium is reachedunder P ∗ = 0 . 01. The other procedure, by which type 2samples are prepared, is meant to represent the oppositesituation, in which aggregation is much faster than com-pression. Thus, while the cell size is fixed and the solidfraction stays equal to Φ I , grains are attributed random(Maxwell-distributed) velocities and left to interact andaggregate until all N of them join to form one uniquecluster. The system is then equilibrated at P ∗ = 0, andcompressed to P ∗ = 0 . 01. To limit the influence of dy-namical effects, the strain rate is requested not to exceeda maximum value ˙ ǫ max during compression. We expressthis condition with the natural inertial time associatedwith the characteristic force F : ( m is the mass of a diskof diameter a ) T = r amF , (4) Sample type No cohesion Type 1 Type 2 N P/K N = 10 − P ∗ = 0 . P ∗ = 0 . 01Φ (no RR) 0 . ± . 001 0 . ± . 001 0 . ± . . ± . 002 0 . ± . 001 0 . ± . defining a dimensionless inertia parameter I a = ˙ ǫ max T . (5) I a is set to 0 . 05 in our simulations. The main set ofsamples of types 1 and 2 (the latter coinciding with “se-ries A” in paper I), to which some non-cohesive onesare added for comparison, is listed in Table II, in whichthe number of available configurations of different sizesis provided, along with solid fraction under the lowestnonzero pressure. All configurations are prepared bothwith ( µ R /a = 0 . µ R /a = 0) RR, withthe parameters of Table I. The initial solid fraction isΦ I = 0 . 36. Type 2 systems are also available under P ∗ = 0, right at the end of the aggregation stage [36],but we regard this intermediate stage as part of the initialpacking process and focus our study on higher pressures(as apparent in Table II, the compression from zero pres-sure to P ∗ = 0 . 01 involves a large density increase, andimportant changes of the microstructure are reported inpaper I). Distant interactions between grain pairs sep-arated by a gap smaller than D are scarce, and “rat-tlers”, i.e. , isolated, free grains with no interactions, areabsent in cohesive systems because of the initial aggre-gation process. Coordination numbers under P ∗ = 0 . z ≃ . z ≃ . V tothe characteristic velocity V ∗ defined in (3). V /V ∗ is setto 9 . V /V ∗ was shown in paper I to have a strong influence onthe initial coordination number z at P ∗ = 0 in sampleswith RR: whereas z is larger than 3 for V /V ∗ = 100, itapproaches 2 for small V , of order V ∗ / 10, in which casethe loopless structures of geometric ballistic aggregationmodels are retrieved. However, this effect is strongly re-duced after the compression step to P ∗ = 0 . D. Simulation procedures 1. Equilibrium conditions One of the specificities of our simulations of cohesivepackings under varying pressure is the approach, com-puting cost permitting, of the quasistatic material re-sponse, in which all configurations remain close to me-chanical equilibrium. Equilibrium conditions have to bestringent enough to enable an unambiguous identifica-tion of the force-carrying contact network and a studyof its elastic properties. Due to the frequent occurrenceof small contact force values, this requires forces to bal-ance with sufficient accuracy. We used similar criteriaas in paper I, which, in agreement with other studies oncohesionless systems [38, 40], were observed to provideadequately accurate force values. The tolerance levels onforce and torque balance equations is expressed in termsof a typical intergranular force value F = max( F , P a ).A configuration is deemed equilibrated when (1) the netforce on each disk is lower than 10 − F ; (2) the total mo-ment on each disk is lower than 10 − F a ; (3) the differ-ence between imposed and measured stresses is less than10 − F /a ; and (4) the kinetic energy per grain is lessthan 5 × − F a . Those conditions being met, we couldcheck that, in the absence of external perturbations (andof thermal motion), no remaining slow motion, creep oraging phenomena were present in our systems: on waitinglonger, only a very slow decrease of the remaining kineticenergy is observed. Furthermore, the computation of thestiffness (or “dynamical”) matrix, see Sec. II D 3 providesan additional stability check. 2. Compression The sample series of Table II are subjected to a step-wise compression cycle. In each compression step, ex-ternal reduced pressure P ∗ is multiplied a constant fac-tor 10 / ≃ . I a = 0 . 05. Parameter I a , on replacing, in its def-inition, F by the force scale aP (in 2D) correspondingto the confining pressure is analogous to inertia param-eter I used to assess dynamical effects in steady shearflow [34, 35], or in the compression of non-cohesive gran-ular packings [38, 41]. The compression program is pur-sued until P ∗ reaches the maximum value 13 . 33, abovewhich negligible plastic collapse is observed. It shouldbe noted that, thanks to the high value of stiffness pa-rameter κ (see Sec. II A), the typical contact deflection aP/K N at this highest pressure level is still very small.Then, the effect of decreasing P ∗ back from its highestvalue to 0 . 01 is also simulated. As no large structuralchanges occur on decompressing the system, larger pres-sure jumps can be imposed on unloading.The simulations are computationally costly, as in somepressure steps equilibration times of order 100 T are re-quired, while the time step for the integration of the equa-tions of motion is a small fraction of p m/K N = T / √ κ .This limits the size and the number of samples, and theuse of small strain rates. Some tests of statistical signifi-cance and rate dependence of the results will be reportedin Section III. 3. Computation of elastic moduli We observe that once samples are equilibrated accord-ing to the conditions of Section II D 1, then the Coulombcriterion | F T | ≤ µF eN , as well as the rolling friction con-dition | Γ R | ≤ µ R F eN are satisfied as strict inequalities inall contacts. No contact is ready to yield in sliding, andwith RR no contact is ready to yield in rolling either.This ensures that the response to small enough exter-nal load increments about a well-equilibrated state willbe elastic and reversible. Elastic moduli express elasticresponse, i.e. , with no effect of tangential or rotationalsliding and no change in contact network topology andgeometry. To compute elastic moduli, we build the stiff-ness matrix K of the contact structure (also taking intoaccount the distant interactions). K [36] is a square ma-trix of order 3 N + 2 (the number of degrees of freedom inthe system), depending on stiffness coefficients K N (re-placed by − F /D for the rare distant attractive bonds), K T , K R (with RR), and on network geometry. K issymmetric, positive definite (once the free translationalmotions of the whole sample as one rigid body are elim-inated) – and thus the stability of equilibrium states ischecked. To compute elastic moduli, one solves a linearsystem of equations: K · U = F ext (6)for the unknown displacement vector U , containing allparticle displacements and rotations, as well as strains( ǫ α ) α =1 , . The right-hand-side of (6) contains externalforces and torques applied to the grains, which are set tozero, and stress increments (∆ σ α ) α =1 , (the same proce-dure is followed in [42] with 2D disk packings and in [43]with 3D sphere packings). On setting ∆ σ = 1, ∆ σ = 0, or vice-versa, one thus gets two separate measurementsof the compliance matrix in our (statistically) isotropicsystems, from which moduli C and C are deduced,and hence the bulk modulus B = ( C + C ) / G = ( C − C ) / III. MATERIAL BEHAVIOR UNDERISOTROPIC LOADA. Compression and pressure cycle withnon-cohesive material Non-cohesive systems of Table II, initially obtainedby isotropic compression of a granular gas (like the 3Dsphere packings of e.g. , Refs. [38] and [44]), are subjectedto a compression cycle, in which reduced pressure P/K N increases from its initial value P /K N = 10 − , up to P /K N = 1 . × − , and decreases back to 10 − .Typical results for the density of systems with andwithout RR are shown on Fig. 1. Changes of solid frac- FIG. 1: (Color online) Φ versus P/K N in pressure cycle with1400 disk samples with and without RR. Blue dashed linescorrespond to elastic response evaluated with the bulk mod-ulus from initial and highest pressure states. tion are very small (of order 10 − , i.e., of order P/K N for the largest pressure), and nearly reversible (morethan 90% of the density increase is recovered on de-compressing), as observed in Ref. [41] with 3D spherepackings. The slight increase of bulk modulus as a func-tion of Φ is due to the larger density of contacts underhigher pressures. One typical feature of frictional, co-hesionless grain packs assembled by direct compressionis the existence of a non-negligible population of “rat-tlers”, i.e., particles that transmit no force (as observede.g. in Ref. [38] in 3D, or Ref. [42] in 2D systems). Thefraction of rattlers x thus exceeds 20% of the grainsunder P in systems with RR in the present case, andreaches 17% without RR. x is reduced to 14% under P/K N = 10 − . The backbone (force-carrying structure)is the set of non-rattler grains, characterized by coordi-nation number z ∗ = z/ (1 − x ) [38]. z ∗ increases with P , as rattlers get captured by the backbone and gapsseparating neighboring grains close in compression.Changes of x and z ∗ are reversed on unloading (withsome moderate hysteresis effect). The increase of z ∗ as afunction of P , above a minimum value z ∗ , which wouldcorrespond to P = 0, is sometimes described by a powerlaw [45]. With such a fit we can estimate z ∗ , and weobtain values close to 3 with RR and about 3 . 12 withoutRR. z ∗ varies by about 10% in the studied pressure in-terval. As in other simulations [38, 46, 47], the minimumcoordination numbers stay above the “critical” value forrigidity, which is equal to 3 without RR and to 2 withRR [36].Cohesionless systems under isotropic pressure cyclesthus behave nearly elastically in an isotropic pressure cy-cle. As the pressure increases by more than 2 orders ofmagnitude, while remaining in the rigid limit of κ ≫ B. Compressing cohesive systems: generalobservations Once subjected to a pressure cycle, as specified inSec II D 2, the material prepared in initially loose states(type 2 of Table II) behaves as shown in Figs 2, 3 and 4.As the pressure increases, so does the density, and thelarge pores present under low P ∗ gradually disappear.The maximum packing fraction, Φ max = 0 . ± . max is smallerthan the solid fraction of cohesionless systems (for whichΦ > . P ∗ ) curves at growing P ∗ , threeregimes can be distinguished. At first, in a range ofreduced pressure P ∗ of the order of the first nonzerovalue (10 − ), thereafter called regime I, Φ remain ap-proximately constant: the contact network supports thegrowing pressure without rearranging. Then, in a sec-ond pressure interval which we shall refer to as regime II,a fast compression is observed. Density variations slowdown in regime III, for P ∗ of order unity, as a maximumsolid fraction Φ max is approached. On reducing the pres-sure, Φ then remains very close to Φ max : the compactionis irreversible.The consolidation curve is similar to the ones obtainedby numerical simulations in Refs. [29, 31], on impos-ing uniaxial strains to loose packings prepared by ananisotropic ballistic aggregation process, although ourstudy differs from these works in several respects (seeSection I). Refs. [29, 31] focus on regime III, and on dy-namical compaction processes, with a shock wave prop-agating through the sample. The variations of solid FIG. 2: (Color online) Equilibrium configuration of a sampleof 1400 disks with RR in initial state, under P ∗ = 0 . 01, forwhich Φ = 0 . . P ∗ = 0 . 178 (different length and forceunits). fraction Φ versus P ∗ are shown in Fig. 5, for three sam-ples of different sizes. Since all three curves are close toone another, we conclude that the macroscopic behav-ior is correctly captured in our simulations. Our resultsfor Φ( P ∗ ) also resemble experimental curves obtained ondifferent materials, such as metallic powders [5], or xero- FIG. 4: (Color online) Same sample as on Figs. 2 and 3,under the maximum pressure P ∗ = 13 . 3. Solid fraction isΦ = 0 . -2 -1 N = 1400 N = 5600 N = 10976 ( % ) P* FIG. 5: Consolidation and decompression curves in 3 samples(with RR) with different numbers of grains, as indicated. graphic toner [10, 48], at least in regimes I and II. Poquil-lon et al. [5], in particular, in an experimental study of ametallic powder, explicitly distinguish three compactionregimes, with the material elastically resisting compres-sion in regime I, and then some plastic compaction, firstattributed to particle rearrangement, as we observe, andlater to contact plasticity. This latter effect, which isnot included in our model, is likely to explain the dif-ference under high P ∗ between many experiments andour results: experimental curves do not appear to ap-proach an asymptotic density, but witness ongoing com-paction up to the highest investigated pressure levels.In the case of metallic powders [5], quite high pressuresare applied (hundreds of MPa), and, as revealed by di-rect microscopic observations, particles fusing or sinter- ing gradually form compact solids. For metal particleswith d=10 µ m diameter, one can estimate the pressure F /d corresponding to P ∗ = 1 to be in the 0 . P ∗ values in the compactionexperiment reveal a different physical origin of densityincrease. The stiffness parameter, κ , is also significantlysmaller in such experiments, with the consequence thatplastic phenomena cannot be ignored (for a definition anddiscussion of κ in Hertzian sphere packings, see [41]).Contact plasticity dominates in the numerical studies ofMartin et al. [18, 19, 20, 21], which focus on very highdensities (beyond the random close packing value), whenthe material, due to sintering, turns into a porous com-pact. Hence only the early stages of metal powder com-paction, in which densities are quite low [5] correspondto our simulations. In the case of the xerographic ton-ers studied in [10, 48, 49], P ∗ = 1, as discussed in [36],rather correspond to P ∼ 10 Pa. Nevertheless, the con-tact behavior, as investigated by atomic force microscopy,is likely to involve plastic effects [50, 51, 52, 53]. C. Regime I: role of the initial assembling process As shown in [36] (paper I), and briefly recalled inSec. II C, assembling conditions have a considerable in-fluence on packing density and microstructure under low P ∗ . It should be assessed to what extent those importantdifferences in the initial configurations affect the plasticconsolidation curve, and whether such a variability tendsto disappear once the material undergoes significant com-paction. This issue is investigated in this section, inwhich the effects of various features of the preparationprocess are observed. The role of some micromechanicalparameters is also discussed. 1. Compaction and aggregation in the assembling stage The most important feature of the assembling processis the competition between compression and aggregation,which leads to the difference between systems of type 1and 2, as defined in [36] and recalled in Section II C. Type1 samples reach a considerably higher densities from thebeginning, under low P ∗ . Fig. 6 compares the subsequentconsolidation curves. As type 1 systems are initially con-siderably denser, they are able to support larger pres-sures before rearranging, hence a wider regime I plateau.However, the pressure increase eventually reaches a highenough value to induce further compaction, and the con-solidation curve is then very close to that of type 2 sys-tems (the difference is actually smaller than the sam-ple to sample r.m.s. fluctuation). Within the accuracyand statistical uncertainty of our simulations, the differ-ence between initial states of types 1 and 2, althoughlarge, thus appears to disappear eventually upon plasti-cally compacting the material. -2 -1 ( % ) P* Type 1, no RR Type 1, RR Type 2, no RR Type 2, RR FIG. 6: Consolidation curve in type 1 and type 2 samples. 2. Effects of first compression step and strain rate In paper I [36] important changes between P ∗ = 0and P ∗ = 0 . 01 in type 2 configurations were reported,as solid fraction Φ increases from Φ I = 0 . 36 to about0.5 (see Table II). One way to limit the effects of thisfirst compression step causes the most dramatic changeis to reduce the strain rate, setting parameter I a to alower value. As shown on Fig. 7, displaying the con-solidation curve obtained in N = 1400 systems with theusual value I a = 0 . 05 and with the smaller one I a = 0 . P ∗ = 0 . 01 is prepared, result ina lower density and tends to turn the initial plateau ofthe Φ( P ∗ ) curve into a gentle ascending slope. Lateron, as consolidation proceeds, very similar curves areobtained with both values of maximum dimensionlessstrain rate I a (Fig. 7), although the smaller error bars(representing sample to sample r.m.s. fluctuations) wit-ness smoother changes and better reproducibility for theslower compression. It may thus be concluded that thequasistatic consolidation curve is quite reasonably ap-proached with the standard compression procedure de-tailed in Section II D 2, for which I a = 0 . 3. Effect of initial agitation and influence of RR The initial agitation velocity (or “granular temper-ature”), as expressed by ratio V /V ∗ in the aggrega-tion stage strongly influences the coordination number.Figs. 8 and 9 show how this initial influence affects thebeginning of consolidation curves and, once again, fadesout later on. Consolidation curves are shown in Fig. 8 fortwo different values of V /V ∗ , one tenfold as large as thestandard value 9 . V on coordination number. An increase ofrolling resistance (with µ R = 0 . . V , stabilizes looser systems under -2 -1 ( % ) P* R / a = 0.005V / V* = 9.5 I a = 0.01 I a = 0.05 FIG. 7: Consolidation curve with two different values of I a . -2 -1 R / a V / V* 0.005 95 0.005 0.095 0.5 95 0.5 0.095 ( % ) P* FIG. 8: Consolidation curve: effect of initial agitation levelin aggregation stage, and influence of RR parameter. low P ∗ , with smaller coordination numbers. However,such a change in material properties does not only affectthe initial, regime I part of the consolidation curve; italso alters the macroscopic mechanical behavior at largerdensities: the slope of the consolidation curve is lower forlarger RR. 4. Conclusion on initial states and regime I Fragile tenuous structures due to aggregation are easilyperturbed and sensitive to many factors in low consoli-dation states. In general, all perturbations favor somekind of preconsolidation effect, inducing denser, bettercoordinated structures. These effects are reduced in eachone of the following situations: (1) if one waits until largeaggregates form before applying a confining pressure; (2)if the initial agitation velocity V is decreased; (3) forslower compression processes, especially when the veryfirst non-vanishing pressure value is imposed; (4) withlarger RR levels. As the material is further compressed -2 -1 R / a V / V* 0.005 95 0.005 0.095 0.5 95 0.5 0.095 z P* FIG. 9: Same as Fig. 8, for coordination number z as a func-tion of P ∗ . in (nearly) quasistatic conditions, the same macroscopicbehavior is retrieved for given microscopic force laws [i.e.,in cases (1) to (3)], irrespective of the initial perturba-tions affecting the beginning of the consolidation process.Though we did not vary the level of viscous dissipationin normal collisions, lower values are expected to inducelarger inertial effects, similarly to a faster compression.On the other hand, viscous forces slowing down the mo-tion of grains relatively to a surrounding fluid (often animportant physical effect in fine powders) could reducethe effects of the initial agitation.Regime I, with no plastic strain, is also observed insome experiments. For example, the response in uniaxialcompression (i.e., σ > σ = σ = 0) of loose aggre-gates of micrometer-sized silica beads assembled by bal-listic deposition – in that case, an anisotropic process inwhich particles are thrown onto a substrate – was studiedby Blum and Schr¨apler [54]. The deposit, with volumefraction Φ ≃ . 15, resists a stress of 500 Pa before plasticcompaction is observed, which corresponds to a “reducedstress”, defined, in analogy with P ∗ , as σ ∗ ≡ σ a /F oforder 10 − . In the simulations of Wolf et al. [31] somefinite initial pressure increment also has to be appliedbefore plastic collapse is observed. D. Regimes II and III:intrinsic consolidation behavior Once the peculiarities of the sample preparation andfirst compression stage are erased, we refer to the mate-rial evolutionas the intrinsic consolidation behavior. Inorder to compare the shape of the consolidation curve toother observations more directly and quantitatively (andalso for a more fundamental reason to be stated further)we subsequently describe it with 1 / Φ, instead of Φ, as afunction of log P ∗ . This conforms to its traditional pre-sentation in the literature [1, 3, 4, 5, 10], which often usesthe void ratio , e = (1 / Φ) − 1. Once the regime I ends, we obtain linear variations of e or 1 / Φ with log P ∗ :1Φ = 1Φ − λ ln P ∗ P ∗ (7)where P ∗ and the corresponding solid fraction Φ are thecoordinates of the point where the system behavior joinsthe intrinsic consolidation curve in the available samples.Parameter λ , known as the plasticity index, is observedin our case to decrease as µ R increases from zero (Fig. 8).We have also observed that the value of this index is notaffected by the friction coefficient: in that sense, µ justdisplaces the whole consolidation curve vertically [53].As the maximum solid fraction Φ max is approached,Eq. (7) is no longer valid, and the asymptotic regime isbetter described with a power law, as in [31]:1Φ = 1Φ max + A ( P ∗ ) α , (8)with a constant A and an exponent α (close to 1 in ourresults). In order to describe the consolidation curve inregimes II and III with a unique functional form, we usethe following relation:1Φ = 1Φ − λ ln ( P ∗ P ∗ (cid:20) − exp (cid:18) − (cid:20) P ∗ P ∗ (cid:21) α (cid:19)(cid:21) /α ) , (9)which introduces additional parameters P ∗ and α , andcrosses over from Eq. (7), for P ∗ ≪ P ∗ , to Eq. (8), for P ∗ ≫ P ∗ . Constant A in (8) is set to λ/ (2 α ) on using (9)for large P ∗ values, and P ∗ is directly related to Φ max :ln P ∗ P ∗ = 1Φ − max . Fig. 10 summarizes the definition and the role of all pa-rameters of relation (9). A fit of our data to relation (9) Region III1/ max ln(P *)ln(P *) ln(P*) Region I Region II Region I: 1/ ~ function(preparation) Region II: 1/ - ln(P*/P* )Region III: 1/ max + const (P*/P* ) - FIG. 10: Schematic view of intrinsic consolidation curve withregimes II and III, and role of parameters introduced inEq. (9). -2 -1 R / a = 0 R / a = 0.005 / P* FIG. 11: (Color online) Consolidation data and fit to Eq. (9),for systems with and without (small) RR. is shown in Fig. 11. It should be noted that even asmall level of rolling resistance changes the plasticity in-dex. Values of parameters are listed in Table III, wherewe also included the fit parameters for the sample with µ R /a = 0 . µ R /a P ∗ Φ λ Φ max α . . 469 0 . ± . 019 0 . . ± . . 005 0 . . 515 0 . ± . 004 0 . . ± . . . . 382 0 . ± . 01 0 . 724 0 . ± . λ , Φ max and α used to fitthe consolidation curve in systems of Table II, and in a samplewith larger RR, with Eq. (9). Correspondingly, P ∗ values are0 . ± . 033 without RR, 0 . ± . 064 for µ R /a = 0 . . ± . µ R /a = 0 . As the consolidation curve in region II, defined by pa-rameters λ and P ∗ , is observed not to depend on initialconditions, our simulations support the following inter-pretation: sooner or later in the process of quasistaticisotropic compression, the system joins, in the P ∗ − Φplane, a certain locus, corresponding to compressive plas-tic yielding. This locus, which acts as an attractor inisotropic compression, is a straight line on using coor-dinates ln P ∗ and 1 / Φ. The value of P ∗ simply signalswhere, depending on the preparation process, the yieldlocus is reached. Table I gives the values of the parame-ters defining the intrinsic curve, and of pressure P ∗ whereit is first reached in type 2 systems of Table II.Consequently, in a system prepared at a lower density,it should be possible to observe a wider interval of the in-trinsic consolidation line. We could explicitly check thisproperty in the case of one sample with N = 5600, forwhich the first nonzero equilibrium confining pressure inthe loading history is equal to 2 × − instead of 10 − .This sample appears to have reached regime II sooner (around P ∗ = 10 − , or possibly below). The correspond-ing data points lie on the intrinsic consolidation curve(or, at least, within a distance smaller than error bars)identified on fitting the data of the main sample series,which had a larger first compression step (to P ∗ = 10 − )and a larger value of P ∗ (about 3 × − ). The yieldlocus can thus be extrapolated to lower pressures anddensities, with the same plasticity index λ . On assem- -3 -2 -1 / P* FIG. 12: Comparison of data obtained on the one low P ∗ sample (open triangles), and Eq. (7) (continuous line) withthe parameters of Table III, as deduced from a fit of the data(black triangles) from the more systematic simulation serieswith larger P ∗ . bling cohesive aggregates with arbitrarily low densities,and on stabilizing them under very low initial pressures,it is conceivable (although increasingly difficult in nu-merical simulation because of the computational cost, aswell as in experiments, because of the system sensitivityto perturbations) to create equilibrium structures withsmaller and smaller densities and to explore an increas-ingly larger interval of the intrinsic consolidation curvein the limit of P ∗ → 0. The corresponding solid fractionΦ would then also tend to zero. This limit is compatiblewith the functional form used in Eq. (7), while the use ofthe alternative form [10, 49],Φ − Φ = ν ln P ∗ P ∗ , would lead to contradictions in the limit of P ∗ → E. Unloading behavior On the Φ versus P ∗ curves we have been showing sofar, that the unloading branch, down to P ∗ = 0 . 01, showsvery little density change. This property is actually sat-isfied on decreasing the pressure from other configura-tions in the compression process. Thus Fig. 13 shows1 -2 -1 Path E ’ DCBA P* ( % ) 1 A B B ’ 2 A C C ’ 3 A D D ’ 4 A E E ’ 5 C ’ C E ’ EB ’ C ’ D ’ FIG. 13: (Color online) Effect of different (isotropic) unload-ing/reloading histories on solid fraction. The direct consoli-dation curve with decompression from the highest pressure,as shown in previous sections, is ABCDEE’ (path 4). Onunloading along lines BB’, CC’, DD’, the system does notrearrange. Such paths are reversible and do not alter the ma-terial state, since paths 4 (small black dots) and 5 (large, openpink circles) superimpose in P ∗ , Φ plane. that, if P ∗ is reduced to the initial level 0 . 01 from differ-ent states on the consolidation curve, density changes arehardly noticeable, and Φ stays very close to the maximumvalue reached at the largest imposed pressure P ∗ c in thepast. Furthermore, it is checked (in the case of sequence4, drawn with open circles in Fig. 13) that the materialmight be reloaded, with no notable density change untilpressure P ∗ c is reached. P ∗ c is known in soil mechanics asthe consolidation pressure , and a material in a state suchthat P ∗ < P ∗ c is said to be overconsolidated . Upon in-creasing the pressure beyond the consolidation value P ∗ c ,the density irreversibly increases, and this compactionis described by the same curve as in the absence of in-termediate pressure cycle: the recompression curve fromC’ retraces back the same evolution from D to E. Thusthe material behavior conforms to the plasticity of claysin isotropic compression [1]. All decompressing paths inthe P ∗ , Φ plane, along which P ∗ < P ∗ c , are reversible.More precisely, they are similar to the pressure cycles ap-plied to cohesionless systems (Fig. 1), and they do notdepart much from the linear elastic response, as shownon Fig. 14.For the largest P ∗ values, adhesion forces are domi-nated by the confining stress and are nearly negligible: onsetting F to zero in equilibrated systems under P ∗ > FIG. 14: (Color online) Analog of Fig. 1, for the unloadingbehavior of a sample with RR from P ∗ = 13 . P ∗ = 0 . P ∗ = 0 . IV. CONSOLIDATION AND DENSITYCORRELATIONS The gradual collapse of the initially open structure ofloose systems, as visually apparent on Figs. 2, 3, and 4and witnessed by the consolidation curve studied in Sec-tion III, can be characterized by the density correlationindicators introduced in paper I.The initial aggregation process was shown in paperI to result in a fractal structure of the density fieldover intermediate scales, between the grain diameter andsome characteristic correlation length ξ . In the presenceof rolling resistance, even with the small value 0 . a adopted for µ R , the observed fractal dimension is com-patible with the result of the ballistic aggregation model, d F ≃ . 55. The ballistic aggregation model is purelygeometric, and corresponds to the irreversible bondingof particles or aggregates in each collision, with contactsthat are rigid in translation and rotation. This limit case,for which the coordination number is equal to 2, is ap-proached under low pressure [36] with large RR or small V /V ∗ . Better coordinated systems obtained with smallRR and/or larger V /V ∗ have the same fractal dimen-sion. Systems with no RR, on the other hand, are closerto dense objets with d F ≃ . ξ is a well-known geometric necessity in alarge system with finite particle packing fraction Φ, be-cause (in 2D) a fractal structure of dimension d F < L exhibits an appar-ent density proportional to L d F − . In physically relevantcircumstances, systems with a finite packing fraction Φand a fractal structure over some distance range have afinite correlation length ξ above which the average value2of Φ is observed. One then has Φ ∝ ξ d F − or ξ ∝ Φ − / (2 − d F ) , (10)the prefactor being specific to the particular system stud-ied. Systems with size L ≫ ξ can then be regarded ashomogeneous packings of fractal “blobs” of (linear) size ξ .Such ideas are quite generally used, and were applied tosemi-dilute polymer solutions [55], to silica [56] or poly-meric [57] gels, in computer simulations of aggregationmodels [58], and to various complex, supramolecular ob-jects like fat crystals [59] or asphaltene aggregates [60].One may expect that the density increase caused by thecollapse, under growing load, of the tenuous structuresformed by cohesive packings corresponds to a decreasein the fractal blob size ξ , while dimension d F still de-scribes the scaling of density correlation at smaller scale.One should then observe the scaling predicted in (10).This implicitly assumes that the small scale structure ofthe packing is not affected by the compaction process,which essentially breaks long, thin junctions and fills thelargest pores. A clue in favor of such a scenario is pro-vided by the results of Sec. III C, which suggest that thesame structure is obtained if the material is directly pre-pared with some value of Φ, or if it is assembled first ina looser state and then isotropically compressed, up tosolid fraction Φ.To compute d F and ξ , we measure the “scattering in-tensity” I ( k ), i.e. the Fourier transform of the densityautocorrelation function, as we briefly recall now (see pa-per I for more details). Density field χ ( r ), taking values1 within particles and 0 outside, is first discretized on aregular mesh, then Fourier transformed, thereby obtain-ing ˆ χ ( k ). We then evaluate I ( k ) = | ˆ χ ( k ) | /A , A beingthe cell surface area. Invoking isotropy, it is a functionof k = || k || alone. I ( k ) should then vary proportionallyto k − d F for a ≪ π/k ≪ ξ , and reach some plateau for k < π/ξ .This approach was used in paper I, and yielded thesame fractal dimension, d F ≃ . 52 in systems with RR,under P ∗ = 0 (solid fraction Φ I = 0 . 36) and Φ = 0 . = 0 . ± . ξ decreased from ξ I = 9 . ± . ξ = 5 . ± . 2. It should be noted thatthese values are roughly compatible with relation (10) (as( ξ I /ξ ) − d F = 1 . ± . / Φ I = 1 . ± . P ∗ = 0 . P ∗ = 0 . P ∗ = 1. These results are averagedover the four largest samples (with RR) of Table II. Inspite of the error bars, I ( k ) exhibits the expected form, itis approximately constant below some crossover wavevec-tor 2 π/ξ which increases with Φ, and then decreases, withslope − d F on a logarithmic plot. Pressure P ∗ = 0 . 178 isthe largest one for which this latter feature is clearly ob-served, and I ( k ) data corresponding to smaller pressuresare intermediate between P ∗ = 0 . 01 and P ∗ = 0 . π/ξ , which values have been estimatedby means of the fit function for I ( k ) presented in paper I. -1.5 -1.0 -0.5 0.0-1.8-1.2-0.60.0 l og (I( k )) log (k/2 ) P*= 0.01 P*= 0.178 P*= 1 d F = . FIG. 15: Scattering intensity per unit area versus wave vector k . Results are averaged over the four largest samples (withRR) of Table II. The curve corresponding to P ∗ = 1 – a flat, low scatter-ing signal – is typical of dense, homogeneous media withno fractal range for density correlations.In view of the small value of ξ reached in the loos-est configurations (those with P ∗ = 0 studied in paperI), relation (10) is difficult to test from density correla-tion data. Another characteristic length scale for densityinhomogeneities, used in paper I, is the (mass) averagedradius of gyration of pores. It may provide an alternativedefinition of a blob size ξ ′ , proportional to ξ . We observed ξ ′ ≃ ξ at P ∗ = 0 . 01, In fact, this equality works wellunder very low consolidations. However, under higherconfining pressures we have observed that the definitionof ξ ′ gives lower values than ξ . Figure 16 is a plot of ξ ′ as a function of pressure.Despite the restricted fractal range, our observationstherefore confirm the validity of the “fractal blob” model,with a constant d F and a correlation length ξ decreas-ing as consolidation proceeds, until a final, homogeneousstructure similar to that of cohesionless packings (albeitsomewhat looser) is obtained. Other values of d F arelikely to be observed with other assembling processes(such as e.g., diffusion-limited cluster aggregation).Values of ξ and d F do not, however, entirely determinethe mechanical properties of the system. The responseof an aggregate to some mechanical perturbation shoulddepend on its connectivity which, as explicitly shown inpaper I, is independent of its fractal dimension (systemswith different µ R and/or prepared with different values of V /V ∗ have the same d F , but very different coordinationnumbers – see also Section III C 3).Results concerning blob sizes in systems without RR,for which d F ≃ . -2 -1 ’ / a P* no RR RR FIG. 16: Average radius of gyration of pores, ξ ′ , versus P ∗ . other types of loose structures. As an example, inthe slow steady state shear flow of a very similar ma-terial simulated in [35] under normal reduced pressure P ∗ = 1 . × − and shear stress σ ≃ . P anisotropicstructures with Φ ≃ . V. PROPERTIES OF EQUILIBRIUM FORCENETWORKSA. Average normal force Formula (2), as explained in Ref. [15] and in paper I,leads to a simple relation between the average normalforce h F N i in equilibrium, pressure P , solid fraction Φand coordination number z : h F N i = π h d i Pz Φ h d i = 7 πaP z Φ . (11)We observed formula (11) (involving the first and sec-ond moments of the diameter distribution) to be accuratein all simulated states despite some approximations in-volved [36]. However, as stressed in paper I, relation (11)fails to estimate the typical contact forces in the networkunder low P ∗ . Those reach values of order F [16]. Nor-mal forces of both signs (as visible on Fig 2 and Fig 3)coexist and, to a large extent, compensate under low P ∗ . B. Coordination numbers In initial low-pressure states, the coordination number, z , as shown on Fig. 17, is nearly equally shared betweenthe contribution z + of compressive bonds and z − of ten-sile bonds. A small population ( z per grain) of contacts carry forces equal to zero (within the numerical toler-ance for force equilibrium). Those contacts, in which thenormal deflection h takes the equilibrium value h for iso-lated pairs [36], tend to be more numerous in the absenceof applied stress, if the aggregation process avoids thebuilding of hyperstatic (overbraced) structures. Theirnumber is quickly reduced once aggregates made under P = 0 are subjected to some external stress and startrearranging.The population of contacts loaded in compression in-creases along the consolidation curve until it dominatesat large P ∗ . Upon unloading, the initial proportion oftensile forces is first retrieved, and z − is eventually, un-der low P ∗ , larger than z + . FIG. 17: Coordination numbers versus P ∗ in compressioncycle (a) without and (b) with RR. Both plots display, fromtop to bottom, z , z + , z − , and z . The error bars (not shown)are about the size of symbols. Arrows close to the curvesindicate the compression branch of the pressure cycle. The total coordination number increases very little inthe pressure cycle. Our observations thus contradictsome statements in the literature [61, 62] relating z toΦ (even in cohesionless systems, Φ and z can vary inde-pendently [38]). In the course of plastic collapse of loose4structures, as the solid fraction increases by more than50%, we observe the number of contacts to increase by5% in systems without RR, and by 12% with RR. Such asmall variation of z in plastic compression contrasts withthe comparatively very fast change of z in the quasielas-tic compression of cohesionless packings, as observed inSection III A, in which z increases by more than 10% forminute density increases.As to the number of distant attractive interactions, i.e.pairs of neighboring grains separated by a gap smallerthan D (contributing to z − ), it is initially very low (typ-ically 10 in a sample of 5600 particles), and then increaseswith P ∗ but remains below 2% of the total number of in-teractions. C. Distribution of forces Normal force distributions are (roughly) symmetricabout zero in initial states under low P ∗ [16], as shownin Fig. 18. Under low P ∗ , tangential forces of order F are also frequently observed [36], and the angle betweenthe total contact force F and the normal unit vector n isnot constrained by the Coulomb condition, which appliesto F + F n rather than to F . This explains the typi-cal patterns of self-balanced contact forces in small grainclusters, where compressive and tensile forces of order F compensate locally, as might be observed on Fig. 2. TheCoulomb condition applying to F , on the other hand, fa-vors alignments and “force chains”. Self-stressed smallclusters form spontaneously when the disks aggregate,except for large RR and/or small V /V ∗ [36].As consolidation proceeds, under growing P ∗ , normalforce distributions develop a wider positive (compressive)side (Fig. 18), while the finite value for F N = − F ischaracteristic of the failure of bonds in traction. Forceseventually scale proportionally to P ∗ at large P ∗ , like incohesionless systems [41, 42], as shown by Fig. 19. When P ∗ reaches values of several unities, the force distributionis similar to that of cohesionless packings, with an addi-tional dwindling population of tensile contacts (Fig. 17).Force distributions in systems with small RR are quitesimilar to those shown in Figs. 18 and 19. D. Forces in dense, overconsolidated states Upon decompressing to low pressure levels, some largercompressive forces ( F N /F reaching 2 or 3) survive andthe distribution is not symmetric (Fig. 18). Such effectsof overconsolidation on contact forces are considerablylarger than in cohesionless granular materials [41]. As inthe case of cohesionless systems [41], we observed thatthe decompression process tends to be affected by dy-namical effects if it is too fast, and the overconsolida-tion effects on force distributions tend to be erased if toomany contacts open in transient stages. The results per-taining to overconsolidated states shown in Figs. 18, 17 FIG. 18: (Color online) Probability distribution function P ( F N ) of static normal force in contacts, versus F N /F , insystems with no RR, for P ∗ = 0 . 01 (black), P ∗ = 0 . 178 (red), P ∗ = 1 (blue), P ∗ = 2 . 37 (green). Distribution widens aspressure increases as indicated by the arrow. P ( F N ) is alsoshown for P ∗ = 0 . 01 for the overconsolidated state (OCS) atthe end of the pressure cycle (pink dashed line).FIG. 19: (Color online) Positive wing of probability distribu-tion function of rescaled normal forces, F N /P ∗ , in systemswith no RR, under P ∗ = 2 . 37 (black crosses), P ∗ = 5 . 62 (redsquare dots), P ∗ = 13 . and 20 were obtained on simply reversing the stepwisecompression program with the parameters indicated inSection II D 2 (i.e. with as many steps in decompressionas in compression).This final force distribution is similar to the one re-ported by Richefeu et al. [16] in simulations of packings ofwet spherical beads, in which cohesion is due to capillaryforces. After assembling the packing under a finite pres-sure and then decompressing to P = 0, these authors ob-serve that the particles tend to form small domains withonly compressive or only tensile forces. Fig. 20 revealsquite similar patterns in overconsolidated states under P ∗ = 0 . 01, with some predominance of the regions undertension, while compressive forces tend to organize more5 FIG. 20: (Color online) Dense overconsolidated state of asample under P ∗ = 0 . 01 (with RR) at the end of the pressurecycle. Color code as on Fig. 2, with distant attractive forces(for which 0 < h < D ) in blue. often in strong force chains. Tensile contacts are morenumerous than compressive ones after the pressure cycle(Fig. 17).To what extent overconsolidation effects on innerstates influence the mechanical properties of cohesivegranular materials (e.g., their response to shear stress)would deserve to be investigated. E. Effect of a large rolling resistance As reported in the previous paragraphs, the small levelof rolling resistance used in most simulations reportedhere ( µ R /a = 0 . λ (see Table III) and changes frac-tal dimension d F (Section IV).In order to understand the mechanisms by which RRaffects macroscopic behavior and geometry, we investi-gated the effects of a large RR ( µ R /a = 0 . 5) in a few1400-disk configurations. Rolling resistance favors forcetransmission along thin strands of particles, each of themin contact with two neighbors (such structures are shownin paper 1 [36, Fig. 20]). Single particle chains are eas-ily disrupted if µ R /a is small, but are quite frequent forsuch RR levels. While single particle chains are easily dis-rupted if µ R /a is small, they become much more frequentfor large rolling resistance Thus coordination numbersmay approach 2 (see Fig. 9). The density and the lengthof such particle chains are also witnessed by the propor- tion x of the contacts that join 2-coordinated disks.Such contacts are impossible in an equilibrium structurewithout RR. x reaches 12% in large RR systems underlow pressure (for Φ in the 0.4 to 0.5 range), down to 1-2% in the main sample series of Table II with small RR( µ R /a = 0 . x . In the limit of z → 2, which isapproached under low pressure for large RR and/or lowvelocity V in the assembling stage, the force network hasa vanishing number of loops and approaches isostaticity,as discussed in paper I [36]. Consequently, as comparedto the case of small or no RR, systems with large RRunder P ∗ ≪ P ∗ of order 10 − , normal forces above F / − F / 10 are extremely scarce (with probability distribu-tion function P ( F N ) in the 10 − range). Furthermore,with P ∗ ∼ 1, while compressive normal forces of order F are frequently observed, P ( F N ) remains below 10 − for F N → − F . This contrasts with the results shown onFig. 18: the proportion of contacts on the verge of ten-sile rupture is much smaller in systems with large rollingresistance. VI. ELASTIC MODULI Elastic moduli are used in experiments [63] and com-puter simulations [43, 64] to express the response of gran-ular materials to small load increments. Their measure-ment, or that of wave velocities, is a non-destructiveprobe of the packing structure. Thus, in the case of co-hesionless bead packings, the simulations of [43] showedthat the moduli are sensitive to coordination number,which can vary independently of the solid fraction, andescapes direct observations [38]. In the present case ofpossibly loose and poorly connected cohesive systems,those moduli also approximately describe the parts ofthe compression curves with no packing rearrangement(Fig. 14), like in cohesionless systems (Fig. 1). A. Elastic moduli of cohesionless packings We first quickly describe the variations of elastic mod-uli in the cohesionless systems of Table II, and their rela-tions to microstructural or micromechanical parameters.Fig. 21 is a plot of bulk and shear moduli versus pressure.Values of moduli are very similar in systems without andwith RR, and vary very slowly with µ R in the latter case.Unlike with Hertzian contacts, local stiffness constants K N , K T do not depend on forces. Consequently, the in-crease of moduli with pressure is moderate. The resultsof Fig. 21 are typical of cohesionless granular systemswith small coordination number [42, 43, 47]. The evolu-tion of bulk modulus is correctly described by the simple6 FIG. 21: (Color online) Bulk and shear moduli of cohesion-less systems (with RR), versus pressure in compression cy-cle. Voigt and Reuss bounds are shown as (red) triangles and(blue) round dots, respectively. Asterisks show values of B obtained on taking a larger rolling stiffness, K R = 10 − K N a instead of K R = 10 − K N a , with the same contact network. estimation formulae recalled below in Sec. VI B, and it isexplained by the increase of coordination number. Shearmodulus G , on the other hand, is somewhat anomalouslylow, witnessing the propensity of a rather poorly con-nected contact network ( z ∗ ≃ . P/K N = 10 − ,without RR) to rearrange under small stress increments, if those are not proportional to the preexisting stresses .The evolution of elastic moduli in the unloading part ofthe pressure cycle (not shown on the figures, for clarity)very nearly reverses the effect of the first compression. B. Simple estimation formulae Bulk and shear moduli are traditionally estimated bythe Voigt or mean field formula [44, 65], which gives up-per bounds [43] B V , G V in terms of contact stiffnessconstants and coordination number z , based on the as-sumption that particle centers move like points of a ho-mogeneously strained continuum. In the present case onehas: B V = z Φ (cid:2) h d i + h d i (cid:3) K N π h d i = 55 z Φ K N πG V = K N + K T K N B V (12)On deriving (12), similar approximations are used asfor (11). The formulae are identical for systems withor without RR, and since we chose K T = K N one hasalso G V = B V .For the bulk modulus, one may also write down a lowerbound B R , the Reuss estimate [43], based on the eval-uation of the elastic energy with trial forces in a load increment. The formula for B R involves moments of thecontact force distribution, specifically the following ratio:˜ Z = h F N + K N K T F T + K N K R Γ ih F N i , (13)in which averages are taken over all contacts carryingstatic normal force F N , tangential force F T and rollingmoment Γ (to be set to zero in the absence of RR). Us-ing (11), one has B R = z Φ h d i K N π h d i ˜ Z = 27 z Φ K N π ˜ Z . (14)This approximation of the bulk modulus becomes exactwhen the force increments caused by an isotropic pressureincrease are proportional to the preexisting forces [43],and hence it tends to be accurate in systems with smalldegrees of force indeterminacy. The ratio of upper tolower bounds for B given by Eqs. (12) and (14) is55 ˜ Z / 54, and the bulk modulus is therefore especiallywell predicted when the force distribution is not toowide [43], and ratio ˜ Z stays close to 1. Thus bulk moduliare rather successfully estimated (see Fig. 21) by B R or B V in the cohesionless case of Sections III A and VI A.Force distributions have often been studied in cohesion-less systems, in which they are strongly constrained bythe no-tension condition, and ˜ Z cannot reach large val-ues ( ˜ Z ≤ . C. Elastic moduli in cohesive packings Elastic moduli as functions of P ∗ during consolidationof cohesive systems are plotted in Fig. 22. Note the log-arithmic scale used for elastic moduli (unlike in Fig. 21).Both bulk and shear moduli are very low at small P ∗ ,which cannot be simply explained by the factor z Φ ap-pearing in estimates (12) and (14) ( z values, see Fig. 17,are similar to those of cohesionless systems while Φ istwice as small at most). Those anomalously low moduliwitness the propensity of the system to rearrange underisotropic as well as under deviatoric stress increments .Moduli in samples with RR (Fig. 22b) have very similarvalues as in the absence of RR, although this may bepartly coincidental, since they are quite sensitive to thevalue of rolling stiffness K R .On decompressing, the moduli (not shown in Fig. 22)stay close to the value reached at the highest pressure.Mean field estimates B V and G V are both too large byfactors of 30 to 50 in loose states. From (11) the averagenormal force h F N i vanishes as P ∗ tends to zero, while thesecond moment is of order F . Moreover, as tangentialforces are not limited by condition | F T | ≤ µF N , but by | F T | ≤ µ ( F N + F ) instead, their contribution to the elas-tic energy is important (and so is that of rolling momentsin systems with RR). Coefficients ˜ Z thus reach values oforder 10 or 10 under low pressure, whence B V /B R ≫ FIG. 22: (Color online) Bulk and shear moduli of cohesivesystems (a) without and (b) with RR, versus (growing) pres-sure. Same symbols and colors as in Fig. 21. which is impossible in cohesionless systems. The Reussbound for B is first (in regime I) too small by a largefactor. Then, it seems to capture the evolution of thebulk modulus in regimes II and III of the consolidationbehavior. Ratio B/B R is reduced to about 2 for P ∗ oforder 0 . 1, and slightly decreases as compression proceeds.It should be recalled, though, that the Reuss formula es-sentially relates the bulk modulus to another unknownquantity, ˜ Z . D. Elastic moduli and force indeterminacy The low value of the shear modulus in poorly coordi-nated cohesionless packings under isotropic stresses (seeFig. 21) has been observed [43, 47] and argued [66] tostem from its tendency to vary proportionally to the de-gree of force indeterminacy per unit area (or volume in3D) when it is small. As the latter (without RR) is pro-portional to ( z ∗ − − x ), one should have G ∗ ≡ G Φ(1 − x ) ∝ z ∗ − . (15) FIG. 23: (Color online) Elastic moduli (no RR) divided by K N Φ(1 − x ), versus z ∗ , the coordination number withoutrattlers. Data with error bars, fitted with the dashed straightline, correspond to G in cohesionless systems. G and B in thecohesive material are respectively shown as (red) crosses andasterisks. Fig. 23 shows our cohesionless packings to abide by thislaw, as the linear variation of G ∗ with z ∗ would predict,within uncertainties, its vanishing for z ∗ = 3. However,it is also obvious from Fig. 23 that the anomalous be-havior of both moduli in loose, cohesive grain assembliesare not simply explained by their low coordination num-ber, except perhaps for the shear moduli of the densestconfigurations (rightmost data points), which, after suffi-cient plastic compaction, become similar to cohesionlesspackings. A coordination number z ∗ just above 3 (with-out RR) characterizes a “barely rigid” contact network,but such a global, average quantity does not account forthe specific heterogeneities of loose cohesive packings. E. Contact forces in a small pressure increment At the microscopic level the elastic response to a smallpressure increment ∆ P determines contact force incre-ments as visualized in Fig. 24. Very strong compres-sive force chains appear, while large parts of the systemcarry very small forces. On sorting the contacts by de-creasing contribution to the elastic energy of the forceincrements balancing ∆ P , less than half of them (46%)contribute 95% of the energy. This proportion increasesto about 65% in the densest configurations, to be com-pared to 68–70% in cohesionless systems. The configu-ration of Fig. 24, in a system with RR, has quite a fewdead ends, i.e. , sets of grains that are connected to therest of the structure but do not belong to any percolatingloop for force (or current) transport through the wholeperiodic cell. With RR, the force-carrying structure coin-cides with the backbone in the sense of ordinary (scalar)percolation theory. The force patterns of Fig. 24 dif-fer from those of Fig. 2, in which the equilibrium forces,prior to the application of ∆ P are shown: some regions,8 FIG. 24: (Color online) Force increments associated with elas-tic response in isotropic compression of system of Fig. 2. Con-tacts are ordered by decreasing contribution to elastic energy,and only the first 46% contact forces corresponding to 95% ofthe energy are drawn (colors as on Figs. 2 to 4). especially the isolated, self-stressed clusters where com-pressions and tensions of order F equilibrate, carry largeforces but are bypassed in the transmission of the pres-sure increment ∆ P . Dead ends contain “islands” of self-balanced forces resulting from the aggregation process,as directly visible on Fig. 2, but they do not participatein the transmission of stress increments and they do notcontribute to elastic moduli.As consolidation proceeds, the repeated application ofpressure increments clearly favors force chains over local-ized self-stressed clusters, and the force pattern adapts tothe external pressure. Hence a closer similarity betweenthe spatial distribution of equilibrium forces under pres-sure P and that of force increments caused by a smallcompression step ∆ P , and a better performance of theReuss estimate. F. Scaling with fractal blob size The inability of the approaches used in cohesionlesssystems to predict the elastic moduli of loose cohesivepackings can be attributed to their ignoring the peculiarnetwork geometry, which is the origin of the strong forceconcentration shown in Fig. 24.In view of the results of Section IV, it is tempting torelate the elastic moduli to the variations of blob size ξ .In scaling arguments about the density, the system canbe regarded as a densely packed assembly of somewhat fuzzy ξ -sized objects, the blobs. To discuss elastic prop-erties, the system is better represented as a network of“superbonds” of length ξ , or effective beams (with whichthe elongated structures carrying stress in Fig. 24 couldbe identified). In such a network, the dominant defor-mation mode is beam bending. The transverse deflection δ in bending of a beam of length ξ , caused by a force F , is proportional to ξ F . Macroscopically, strains areof order ǫ = δ/ξ , while F corresponds to stress σ by F ∝ σξ in 2D. Consequently, the scaling of elastic mod-uli σ/ǫ with length ξ should involve a factor ξ − . (Somepossible corrections to exponent 3 are possible, althoughthe appropriate value in, e.g., the case of percolation net-works of beams is very close to 3 [67, 68]). As ξ varies bya factor of 3 or 4 within the scaling range (see Fig. 16),relation B ∝ ξ − would predict an increase of moduli bya factor of a few tens.Although this can be regarded as a fair estimate (seeFig. 22), it should be admitted that the fractal range isvery likely too restricted for such scaling laws to applywithout important corrections. With sufficiently largerolling resistance, the “beams” can be reduced to singleparticle chains, which, as we now show, enables simpleranalyses of their bending stiffness. G. The case of a large rolling resistance With large RR the prevalence of particle strands asforce-transmitting structures (Sec. V E) influences elas-tic properties. As noted above, linear structures tend todeform like bending beams, with a compliance propor-tional to the third power of their length. In the case ofsingle linear strands, connections with contact propertiesare easily made more explicit. Consider, e.g., a straight,linear chain of n identical disks of radius R , with n − K N , K T and K R . Then, in the elastic regime, all intermediate diskscan be suppressed and the interaction between the ex-treme ones, numbers 1 and n along the chain, can bereplaced by an effective one between two disks of radius( n − R , and compliances 1 /K ( n ) N , 1 /K ( n ) T , and 1 /K ( n ) R for normal, tangential and rolling relative motion, with: K ( n ) N = n − K N K ( n ) T = n − K T + ( n − n − n + 6) R K R K ( n ) R = n − K R (16)For large n the tangential compliance is much larger thanthe longitudinal and rolling ones, so that long chains be-have as beams, which essentially deform in bending. Thelocal bending stiffness EI of the beam (i.e. the productof the material Young modulus by the moment of inertiaof the beam section) corresponding to the chain of parti-9cles in the continuous limit is EI = 2 RK R . (This coeffi-cient expresses the proportionality of bending moment torotation angle gradient). For n ≫ 1, the bending springconstant 3 EI/l (expressing the transverse force to trans-verse deflection relationship) is correctly identified from K ( n ) T given in (16), using the length l = 2( n − R of thestraight n -particle strand.Remarkably, the bending elasticity of small linearstrands of micrometer-sized colloidal particles bound byadhesive forces has recently been measured by means ofoptical tweezers [69]. Colloidal gels of polymer parti-cles [70, 71, 72] should thus be modeled as cohesive par-ticle assemblies with rather large RR level.It is easy to check (consider e.g., two such chains join-ing at their ends at some angle) that for all strand shapesother than straight lines, the extremities will be coupledby spring constants of order K ( n ) T for both longitudinal(parallel to end-to-end vector) and transverse relative dis-placements. Consequently, the macroscopic elastic mod-uli should be proportional to rolling stiffness constant K R . Fig. 25 shows that this proportionality is approx-imately satisfied in the loosest states of a system withrolling friction µ R /a = 0 . 05, in which three different val-ues of K R were used to evaluate the elastic response.Elastic moduli of denser states, however, depart from FIG. 25: (Color online) Bulk (filled symbols) and shear (opensymbols) moduli, normalized by K R , in low pressure statesof a sample for which µ R /a = 0 . 05 and K R = 10 − K N a (black squares). Results obtained on evaluating moduli with K R = 10 − K N a and with K R = 10 − K N a are respectivelyshown as red triangles and blue circles. this behavior. Therefore, the scaling of elastic moduliwith typical strand length (as suggested in Section VI F)is limited to low consolidation states. With small or van-ishing RR, single particle strands are replaced by thickerjunctions, which further restricts the consolidation pres-sure range for which elasticity is dominated by beambending. VII. PLASTIC CONSOLIDATIONMECHANISM: QUALITATIVE ASPECTS Cohesionless granular assemblies, if subjected to stressincrements that are not proportional to initial stresses,essentially deform because the contact network gets re-peatedly broken and repaired [40, 73]. Macroscopicstrains, once they exceed the very small scales associatedwith the response of given contact networks [40, 43], thusresult from a sequence of rearrangement events or micro-scopic instabilities, during which the granular packingloses its coherence and gains some finite amount of kineticenergy, even for arbitrarily slow applied stress changes.Collisions and appearance of new contacts stabilize thepacking at the end of each microscopic rearranging event.This process gradually changes the topology of the con-tact network, and produces specific evolution of its fabric(orientation anisotropy).The mechanism of plastic collapse in isotropic compres-sion of loose cohesive assemblies with small or vanishingRR in contacts, as observed in the present study, is sim-ilar. Just like in cohesionless systems under shear [40],we expect the frequency of occurrence of rearrangements,along the loading path, to increase, and the correspond-ing strain jumps to decrease, as the size of samples grows,and thus the consolidation curve should be smooth inthe thermodynamic limit. Due to the specific geome-try of loose systems, in which dense zones are weaklyconnected through thin arms, better connected, solid-like regions tend to move like rigid bodies, while frag-ile junctions break and rearrange, so that initially largeholes gradually fill up. Fig. 26 illustrates this scenario.Displacements are depicted as arrows, pointing from thecurrent positions to the ones reached in the next equilib-rium configuration in the stepwise compression sequence.The more densely packed, nearly rigid regions (markedwith dotted lines) are easily identified by direct visualinspection. Fig. 26 also shows that the contact networkundergoes relatively small topological changes, as morethan 90% of contacts are conserved. The rate of con-tact change, and the evolution of coordination numberwith strain are significantly smaller than in cohesionlesssystems undergoing, e.g., shear deformation. During thecompaction of loose samples the dense regions collide andslide past one another, along thin sheared zones wheremost of the broken contacts are found.In the case of large RR, the peculiar microstructureinvolving single particle chains might lead to a differ-ent deformation mechanism. Unlike multiply connectedjunctions, simple strands can yield in bending withoutbreaking: they fold at some contact, where the rollingfriction threshold is reached, thereby releasing bendingelasticity. This mechanism is observed in experimentson single chains of colloidal particles [69, 71]. One thusexpects fewer contact losses in plastic compression.To follow more closely the rearrangement sequencesin the course of compaction, it is appropriate to moni-tor changes in the list of contacts during the motion be-0 FIG. 26: Equilibrium particle positions in 1400 disk samplewith small RR under P ∗ = 0 . P ∗ = 0 . 042 are shownas arrows (global density change ∆Φ = 0 . tween two equilibrium configurations. As an example, letus consider the evolution between equilibrated states as P ∗ increases from 0 . 177 to 0 . µ R /a = 0 . µ R /a = 0 . 5) RR. Table IV gives the changes insolid fraction and coordination number, and numbers ofmaintained, destroyed and created contacts in this com-pression step. Successive configurations separated by afixed time interval ∆ t = 0 . T are compared and Fig. 27plots the number of destroyed and created contacts asfunctions of time. For the same strain increment, con-tact losses, as a function of global strain, are significantlyless frequent in the sample with large RR. This fact is re-flected both in the data of Table IV, where global changesare recorded, between the initial and final states, and inthose of Fig. 27, where successive changes over time inter-vals ∆ t are detailed. As a consequence, while the coor-dination number hardly changes during consolidation insystems with small or vanishing RR (see Fig. 17), it grad-ually increases from an initial value close to 2 to nearly 3in systems with large RR (Fig. 9). The lesser importanceof tensile contact rupture in the plastic compression of as-semblies with large RR is also witnessed by the normalforce distribution (Section V E): forces approaching − F are quite scarce, as opposed to the situation in samples µ R /a ∆Φ(%) ∆ z (%) N (=) N ( − ) N (+) z ), and numbers of maintained ( N (=) ),destroyed ( N ( − ) ) and created ( N (+) ) contacts in a 1400 diskssample, with small or large RR, in the compression step be-tween P ∗ = 0 . 177 and P ∗ = 0 . µ R /a = 0 . 005 the proportions x + and x − of gained and of lostcontacts with respect to the previous recorded list are respec-tively shown with red square dots and triangles – the latterbeing connected with a dashed line. A similar code is used for x + and x − values in a sample with large RR ( µ R /a = 0 . without RR (Fig. 18). With small RR, some single par-ticle chains are also present, although shorter and lessnumerous. The sensitivity of plasticity index λ to therolling friction is likely to be explained by different rup-ture mechanisms, the importance of folding rearrange-ments growing with the level of rolling resistance. VIII. CONCLUSION To summarize, we have used numerical simula-tions to observe and characterize, at the macroscopicand microstructural levels, the consolidation behavior,in isotropic compression, of model cohesive powders.Macroscopic constitutives laws for quasistatic loading,unloading and elastic responses were shown to be reason-ably well approached. The material behavior was inves-tigated for a range of densities that is wider than in mostsimulation studies of cohesive granular materials. Theconsolidation process goes through three stages. In a first1regime, which is sensitive to the assembling procedure, noplastic collapse occurs, as the agitation in the assemblingprocess has stabilized a strong enough microstructure towithstand a finite pressure increase. The normal forcedistribution widens until a significant fraction of contactsare on the verge of tensile rupture. The initial system ge-ometry, which changes very little in regime I, is that of adense assembly of fractal blobs, with dimension d F takingthe universal value associated with the aggregation pro-cess (here, ballistic) implemented in the sample prepara-tion stage. The blob size ξ (at most, between 5 and 10grain diameters in the present case) can be identified onstudying density correlations. The subsequent consolida-tion behavior is remarkably independent on initial con-ditions, which merely determine where the intrinsic con-solidation curve in the Φ- P ∗ plane is first met. The samecurve is then followed whatever the initial conditions asthe material is further compressed. This behavior cor-responds, at the microscopic level, to a gradual changeof the blob size. The curve in regime II has the sameshape as reported in the soil mechanics literature, andthe consolidation pressure is a plastic threshold belowwhich the material response is approximately elastic (likethe behavior of a cohesionless granular material underisotropic load). Elastic moduli increase rapidly with con-solidation pressure or density. The cases of small RR orwithout RR should be distinguished from the situation ofstrong rolling resistance, although, in both cases, the mi-crostructure of loose packings might be viewed as denser,better connected regions joined by thin arms. In the firstcase, loose packings collapse when the tensile strengthof contacts is overcome by the externally imposed forces,preferentially within the fragile junctions between adja-cent denser blobs. Systems with strong RR, on the otherhand, contain single particle strands, which tend to foldwithout breaking in plastic compaction. While small RRsystems gain very few contacts in the consolidation pro-cess, the coordination number might increase from nearly2 to 3 with large RR. Eventually, the material approachesa limiting, maximum density (regime III), as the pack-ing structure resembles that of a cohesionless system, for P ∗ ≫ P ∗ is due to plasticdeformation of contacts.The fractal blob size ξ , depending on solid fraction Φ,is a central microstructural feature, based on which somescaling laws for elastic properties can be attempted. Itis also tempting, beyond the qualitative description ofthe microstructural changes associated with the consoli-dation process, to try to predict the consolidation curvefrom such geometric data. Yet scaling laws only apply toa restricted part of the consolidation pressure interval.Our results, in many respects, emphasize importantqualitative differences between cohesive and cohesionlessgranular assemblies. The existence of stable loose struc-tures and the consolidation phenomenon are the mostimportant differences in macroscopic behavior brought about by cohesion. At the microstructural level, unlikein cohesionless packings, the typical values of intergran-ular forces, or the force distribution, are not as simplyestimated in cohesive systems, in which attractive andrepulsive contact forces of the order of tensile strength F tend to compensate under low pressure. In particu-lar, compression cycles stabilize self-balanced force net-works with large compression forces. Unlike in granu-lar packings devoid of cohesion, the coordination numberdoes not appear to be a significant state variable in co-hesive systems with low RR, as it hardly changes alongthe consolidation curve. With large rolling resistance, itwitnesses, however, the formation of loops under com-pression. While cohesionless assemblies with low coordi-nation number usually contain many rattlers, all parti-cles in cohesive packings are connected to the same con-tact structure, which is rigid, but comprises lots of “deadends” or “side arms”, which might bear self-balancedforces but do not participate in the transmission of exter-nal stresses. Some of these new features can be summedup on remarking that loose powders are similar to gelsas much as to granular packings with no cohesion.Our investigations should be pursued in several direc-tions. On the theoretical side, the connections betweenmacroscopic properties and microstructure could be stud-ied more quantitatively. The behavior of loose cohesivepackings under general stress states should be investi-gated. Thus one may determine whether such constitu-tive laws as the Cam-clay model [1] apply to the simu-lated material. And finally, more quantitative agreementwith experiments and real materials should be sought.In spite of some obvious steps (e.g., one should simulate3D systems), this latter objective looks daunting. Onemajor difficulty is the importance of hydrodynamic ef-fects at the assembling stage, when the microstructureand the fractal dimension of aggregates are determined.While we have bypassed this problem on implementingballistic aggregation, it is necessary to investigate thebehavior other possible kinds of aggregates, by dealingwith some tractable model for hydrodynamic forces. Itis hopefully possible to introduce some mechanics and in-tergranular interactions within the models used with ge-ometric aggregation rules (such as, e.g., diffusion-limitedcluster-cluster aggregation). Then, another difficulty isthat many parameters associated with the contact law(such as friction coefficient, rolling friction, rolling stiff-ness constant) should be identified for a real material tobe investigated at the grain level. In this respect, therecent progress of experimental methods of microscopicinvestigation seems quite promising, as formerly inacces-sible parameters ruling interparticle contact mechanicsare now beginning to be measured in model materials,thanks to particle-scale observation and micromanipula-tion techniques [50, 51, 69, 71, 74].2 Acknowledgments This work has been supported by the Ministerio deEducaci´on y Ciencia of the Spanish Government under project FIS2006-03645 and by the Junta de Andaluc´ıaunder project FQM-421. [1] D. M. Wood, Soil Behaviour and Critical State Soil Me-chanics. (Cambridge University Press, England, 1990).[2] J. K. Mitchell, Fundamentals of soil behavior (Wiley,New York, 1993).[3] J. Biarez and P.-Y. Hicher, Elementary Mechanics of SoilBehaviour (A. A. Balkema, Rotterdam, 1993).[4] J. H. Atkinson, An introduction to mechanics of soilsand foundations: Through critical state soil mechanics. (McGraw-Hill, International Series in Civil Engineering,New York, 1993).[5] D. Poquillon, J. Lemaitre, V. Baco-Carles, P. Tailhades,and J. Lacaze, Powder Tech. , 65 (2002).[6] A. R. Cooper Jr. and L. E. Eaton, Journal of the Amer-ican Ceramic Society , 97 (1962).[7] J. S. Reed, ed., Principles of Ceramic Processing (Wiley,New York, 1995).[8] D. Falgon, E. Vidal-Sall´e, J.-C. Boyer, R. Peczalski, andJ. Andrieu, Powder Technology , 183 (2005).[9] J. G. Mallol Gash, Ph.D. thesis, Institut de Tecnolo-gia Cer`amica (ITC), Universitat Jaume I, Castell´o de laPlana, Spain (2006).[10] A. Castellanos, Advances in Physics , 263 (2005).[11] P. Pierrat and H. S. Caram, Powder Technology , 83(1997).[12] T. Gr¨oger, U. T¨uz¨un, and D. M. Heyes, Powder Technol-ogy , 203 (2003).[13] Z. Fournier, D. Geromichalos, S. Herminhaus, M. M. Ko-honen, F. Mugele, M. Scheel, B. Schulz, C. Schier, R. See-mann, and A. Skudelny, Journal of Physics CondensedMatter , 5477 (2005).[14] P. A. Cundall and O. D. L. Strack, G´eotechnique , 47(1979).[15] V. Richefeu, M. S. El Youssoufi, and F. Radjai, PhysicalReview E (2006).[16] V. Richefeu, F. Radjai, and M. S. El Youssoufi, Euro.Phys. J. E , 359 (2006).[17] C. M. Kong and J. J. Lannutti, Journal of the AmericanCeramic Society , 685690 (2000).[18] C. L. Martin and D. Bouvard, Acta Materialia , 373(2003).[19] C. L. Martin, Acta Materialia , 4589 (2003).[20] C. L. Martin, D. Bouvard, and S. Shimai, Journal of theMechanics and Physics of Solids , 667 (2003).[21] C. L. Martin, Journal of the Mechanics and Physics ofSolids , 1991 (2004).[22] J. Nam and J. J. Lannutti, Journal of the American Ce-ramic Society , 557 (2004).[23] C. L. Martin and D. Bouvard, Journal of the AmericanCeramic Society , 3379 (2006).[24] S. Luding, Powder Technology , 45 (2005).[25] R. Tykhoniuk, J. Tomas, S. Luding, M. Kappl, L. Heim,and H.-J. Butt, Chemical Engineering Science , 504(2007).[26] K. Dong, R. Yang, R. Zou, and A. Yu, Physical ReviewLetters , 145505 (2006). [27] A. Misra, Materials and manufacturing processes , 925(1996).[28] D. Kadau, G. Bartels, L. Brendel, and D. E. Wolf, Com-putational Physics Communications , 190 (2002).[29] D. Kadau, G. Bartels, L. Brendel, and D. E. Wolf, PhaseTrans. , 315 (2003).[30] G. Bartels, T. Unger, D. Kadau, D. E. Wolf, andJ. Kert´esz, Granular Matter , 139 (2005).[31] D. E. Wolf, T. Unger, D. Kadau, and L. Brendel, in Pow-ders and Grains 2005 , edited by R. Garc´ıa Rojo, H. J.Herrmann, and S. McNamara (Balkema, Lisse, 2005), pp.525–533.[32] R. Brewster, G. S. Grest, J. W. Landry, and A. J. Levine,Physical Review E , 061301 (2005).[33] L. Aarons and S. Sundaresan, Powder Technology ,10 (2006).[34] P. Rognon, J.-N. Roux, D. Wolf, M. Naa¨ım, andF. Chevoir, Europhysics Letters , 644 (2006).[35] P. G. Rognon, J.-N. Roux, M. Naa¨ım, and F. Chevoir, J.Fluid Mech. , 21 (2008).[36] F. A. Gilabert, J.-N. Roux, and A. Castellanos, PhysicalReview E , 011303 (2007).[37] M. Parrinello and A. Rahman, Journal of AppliedPhysics , 7182 (1981).[38] I. Agnolin and J.-N. Roux, Physical Review E , 061302(2007).[39] A. Tordesillas and D. C. Stuart, Powder Technology ,106 (2002).[40] J.-N. Roux and G. Combe, C. R. Acad´emie des Sciences(Physique) , 131 (2002).[41] I. Agnolin and J.-N. Roux, Physical Review E , 061303(2007).[42] E. Somfai, J.-N. Roux, J. Snoeijer, M. van Hecke, andW. van Saarloos, Physical Review E , 021301 (2005).[43] I. Agnolin and J.-N. Roux, Physical Review E , 061304(2007).[44] H. A. Makse, N. Gland, D. L. Johnson, and L. M.Schwartz, Phys. Rev. Lett. , 5070 (1999).[45] C. S. O’Hern, L. E. Silbert, A. J. Liu, and S. R. Nagel,Phys. Rev. E , 011306 (2003).[46] H. P. Zhang and H. A. Makse, Physical Review E ,011301 (2005).[47] E. Somfai, M. van Hecke, W. G. Ellenbroek,K. Shundyak, and W. van Saarloos, Physical Review E , 020301 (2007).[48] M. A. S. Quintanilla, Ph.D. thesis, Facultad de F´ısica,Universidad de Sevilla, Seville, Spain (2003).[49] A. Castellanos, J. M. Valverde, and M. A. S. Quintanilla,Physical Review Letters , 075501 (2005).[50] D. M. Schaefer, M. Carpenter, B. Gady, R. Reifenberger,L. P. DeMejo, and D. S. Rimai, J. Adhes. Sci. Technol. , 197 (1994).[51] M. Reitsma, V. Craig, and S. Biggs, Int. J. Adhes. Adhes. , 445 (2000).[52] M. A. S. Quintanilla, A. Castellanos, and J. M. Valverde, Physical Review E , 031301 (2001).[53] F. A. Gilabert, Ph.D. thesis, Facultad de F´ısica, Univer-sidad de Sevilla, Seville, Spain (2007).[54] J. Blum and R. Schr¨apler, Physical Review Letters ,115503 (2004).[55] P.-G. de Gennes, Scaling Concepts in Polymer Physics (Cornell University Press, Ithaca, 1979).[56] G. Dietler, C. Aubert, D. Cannell, and P. Wiltzius, Phys-ical Review Letters , 3117 (1986).[57] J.-P. Cohen-Addad, ed., Physical properties of polymericgels (Wiley, New York, 1996).[58] P. Meakin, Journal of Sol-Gel Science and Technology , 97 (1999).[59] S. S. Narine and A. G. Marangoni, Physical Review E , 1908 (1999).[60] J.-N. Roux, D. Broseta, and B. Dem´e, Langmuir , 5085(2001).[61] E. Jaraiz, S. Kimura, and O. Levenspiel, Powder Tech. , 23 (1992).[62] R. Y. Yang, R. P. Zou, and A. B. Yu, J. Appl. Phys. ,3025 (2003). [63] R. Kuwano and R. J. Jardine, G´eotechnique , 727(2002).[64] H. A. Makse, N. Gland, D. L. Johnson, and L. Schwartz,Physical Review E , 061302 (2004).[65] K. Walton, Journal of Mechanics and Physics of Solids , 213 (1987).[66] M. Wyart, Annales de Physique Fr. , 1 (2006).[67] A. Kantor and I. Webman, Physical Review Letters ,1891 (1984).[68] S. Roux, Journal of Physics A , L351 (1986).[69] J. P. Pantina and E. M. Furst, Physical Review Letters , 138301 (2005).[70] J. P. Pantina and E. M. Furst, Langmuir , 3940 (2004).[71] E. M. Furst and J. P. Pantina, Physical Review E ,050402(R) (2007).[72] J. P. Pantina and E. M. Furst, Langmuir , 1141 (2008).[73] L. Staron, J.-P. Vilotte, and F. Radjai, Physical ReviewLetters , 204302 (2002).[74] L.-O. Heim, H.-J. Butt, J. Blum, and R. Schr¨apler, Gran-ular Matter10