Scaling laws for the mechanics of loose and cohesive granular materials based on Baxter ’ s sticky hard spheres

We have conducted discrete element simulations (PFC3D) of very loose, cohesive, granular assemblies with initial configurations which are drawn from Baxter’s sticky hard sphere (SHS) ensemble. The SHS model is employed as a promising auxiliary means to independently control the coordination number zc of cohesive contacts and particle volume fraction φ of the initial states. We focus on discerning the role of zc and φ for the elastic modulus, failure strength, and the plastic consolidation line under quasistatic, uniaxial compression. We find scaling behavior of the modulus and the strength, which both scale with the cohesive contact density νc = zcφ of the initial state according to a power law. In contrast, the behavior of the plastic consolidation curve is shown to be independent of the initial conditions. Our results show the primary control of the initial contact density on the mechanics of cohesive granular materials for small deformations, which can be conveniently, but not exclusively explored within the SHS-based assembling procedure.


I. INTRODUCTION
Identifying microstructural controls of macroscopic material behavior is key for the understanding, upscaling, and modeling of heterogenous or disordered materials.The subclass of granular materials is often studied by discrete element simulations to address material behavior for different applications.A large amount of work has been done in the past to characterize the mechanics of dense cohesionless materials (e.g., Refs.[1][2][3][4][5][6][7]), mainly in the vicinity of the jamming transition (or the critical state [8]).In contrast, in loose and cohesive granular systems (e.g., [9][10][11][12][13][14]), the mechanical behavior originates from properties of the tenuous, microstructural grain networks that sustain external loads.From the analysis of consolidation, where plastic effects are negligible, it was concluded [11,12] that the mechanical behavior is predominantly influenced by the microstructure inherited from the assembling procedure.
The most prominent microstructural properties of granular systems are particle volume fraction φ and coordination number z c .For dense, noncohesive systems close to jamming transition, both quantities are interrelated, and universality dictates the coordination number to be a function of packing fraction (cf., e.g., [4,15,16]).Introducing cohesive bonds into the system allows generation of assemblies at small volume fraction where the influence of both parameters may decouple [17].Agnolin and Roux highlighted in [18] that a compression-decompression cycle, although almost reversible in terms of density, can significantly lower the coordination number.However, elastic moduli are strongly sensitive to the coordination number [12].These authors suggested that elastic moduli should be rather related to the contact density as an indicator of the internal state.Discerning the role of φ and z c becomes particularly important for loose states, where the number of possible configurations for a given packing fraction increases.
It is thus desirable to augment the control of initial states used in discrete element modeling.Initial sphere configurations are often generated from ballistic deposition [13] or from aggregation schemes based on equilibrium states of hard sphere fluids [12].The exact influence of volume fraction φ and coordination number z c is, however, difficult to assess, partly due to the lack of analytical approximations available for rule-based definitions of the assembling procedures for initial states.
Along these lines it has been previously suggested [19] that Baxter's model of sticky hard spheres (SHS) [20] is an interesting candidate to address generic questions of disordered materials from a particle-based perspective.In the original work [20] SHS were introduced as a minimal particle model for a liquid-gas transition.The SHS model is also a common starting point for other systems, such as short-range attractive, colloidal suspensions [21].Due to the relevance of gelation and glassy phenomena in colloidal systems, the SHS system was early investigated from the viewpoint of topology and connectivity [22].Such a geometric-configurational perspective of the underlying (thermodynamic) equilibrium states reveals that the SHS model undergoes a percolation transition, where mechanically stable configurations can be achieved also at very low volume fraction.The competition between percolation on one hand and phase separation on the other hand [23] gives rise to rich physical behavior with a strong impact on rheology and mechanics.If adopted from a granular perspective, SHS facilitates a new and interesting view on static and dynamic properties of loose cohesive granular materials.A unique feature of Baxter's SHS is the possibility to parametrically sweep through completely different microstructures which have identical packing fraction but different coordination numbers.The appeal of investigating SHS as initial states for cohesive granular assemblies is also enhanced by the analytical progress which has been made to evaluate thermodynamic properties [24].At the level of the Percus-Yevick (PY) approximation, configurational properties of the initial state, such as pair correlations or percolation loci, are known exactly.This is useful despite well-known limitations of PY in predicting the phase diagram of an SHS fluid [23].As mentioned in [25], the PY approximation does "remarkably well" for the coordination number as the key microstructural descriptor of granular systems.
It is the aim of the present paper to exploit the properties of SHS assemblies to discuss the influence of volume fraction φ and coordination number z c on the elasticity and strength of loose and cohesive granular materials.We use discrete element simulations (PFC3D, Itasca [26]) and focus on the quasistatic (uniaxial) compression in the absence of gravity.The elastic modulus and failure strength are shown to be universal functions of the initial cohesive contact density, which we define by ν c := z c φ following the nomenclature in [27,28].On the other hand, plastic effects are shown to be independent of the initial conditions.We believe that our simulations help to pinpoint the unifying aspects of mechanics, failure, and rheology of heterogeneous, foamlike, granular, and soft matter [29,30].

A. Sticky hard spheres
The SHS model is defined by a classical fluid of monodisperse hard spheres with diameter d interacting via a pair potential The model is defined in the limit δ → 0 in which the finite range attraction tends to a contact adhesion controlled by the stickiness parameter τ .In this athermal limit, τ effectively takes the meaning of the temperature such that high values of τ correspond to low adhesion strengths, and nonsticky, hard spheres must be recovered by τ → ∞.

B. Microstructure and coordination number
Thermodynamic and configurational properties of SHS can be investigated within the Percus-Yevick approximation [24], giving rise to closed form expressions in terms of φ and τ .Despite well known shortcomings [23], we are guided by PY to illustrate those SHS microstructures which can be used as meaningful initial states for a cohesive granular material.
To motivate the choice of simulation parameters, we show an overview of the PY phase diagram in the φ − τ plane in Fig. 1 together with the coexistence line and the critical point of the underlying liquid-gas transition.The SHS system exhibits rich microstructural variability above the coexistence line τ coex (φ), where the attractive forces lead to a percolation transition in the high-temperature region.The percolation transition is reflected by the behavior of the average coordination number of SHS, which is well described by PY [25].For a given packing fraction φ and adhesion τ , the PY coordination number z PY c is given by where λ is the smallest solution of the quadratic equation A contour plot of the average coordination number is added in Fig. 1.This indicates that the locus of the percolation points τ perc (φ) in the φ − τ plane, originally identified with the divergence of the mean cluster size [22], corresponds to z PY c = 2, as a mean-field condition of minimum connectivity.
In the region between τ perc (φ) and τ coex (φ) the system is in a gel state in which the competition between percolation and phase separation dictates the configurational properties of SHS [23].In this region the system exhibits both topological clustering by increasing connectedness (percolation transition) and Euclidean clustering by increasing density fluctuations (first-order transition).The existence of a macroscopic cluster justifies the interpretation of the system as a cohesive solid, if (cemented) bonds are introduced between particles in contact.For the discrete element method (DEM) simulations below, we chose the black dots in Fig. 1 to cover a significant variation in model parameters φ and τ .We stress that the original, singularadhesive pair interactions of SHS from Eq. (1) only serve as an auxiliary means to generate initial sphere configurations in which new pair interactions are introduced for the DEM simulations.

C. Monte Carlo simulations of SHS
The implementation in [31] is used to generate the initial configurations of the SHS ensemble.The Monte Carlo simulation algorithm in [31] is a generalization of the algorithm of [ 32,33] to multisize sticky spheres.In this implementation, the Metropolis shuffling scheme [34] is used to create an initial realization of randomly distributed spheres.This initial configuration then undergoes a number of passes to form and break bonds among particles to reach an SHS equilibrium state.Only transitions among four binding states of unbonded (α = 0), single bond (α = 1), double bond (α = 2), and triple bond (α = 3) are considered.The transitions among states are driven by transition probabilities, which are proportional to the effective volumes V α eff associated with the binding state α.The effective volume is assigned to each state following the Kranendonk-Frenkel algorithm [33], as a measure of all the degrees of freedom of a test particle k after making a α bond, which is inversely related to the stickiness parameter τ through a multiplicative factor of (d/12τ ) α , where d is the diameter of the sphere.The total effective volume for the binding state α of a particle k is simply the summation of all the partial effective volumes of all possible combinations of particle k and other particles that can form a α bond.The trial transition of particles is carried for each sphere in the ensemble one by one in every Monte Carlo simulation cycle.If the particle is currently in contact with more than three spheres or its transition location to be determined turns out overlapping with one of the other spheres in the system, no displacement will be made.To find the transition location, we first choose a binding state for the test particle, disregarding what current binding state it is in.We then select a specific particle combination to form the binding state other than being unbonded.This can be a candidate particle, or pair, or triplets to form a single, double, or triple bond, respectively, and the final step is to determine an exact location to place the test particle.This can be done by choosing the position of the test particle center on the spherical surface of radius d centered at the candidate particle to form a single bond, or selecting the center location of the test particle on the circle that can form a double bond with the candidate pair, or to choose between above or under the surface defined by the centers of the candidate triplets in the case of a triple bond.All possible choices are equally treated in these three steps.If the next state is chosen to be unbonded, the particle is then equal-probably transited to any location in the simulation volume subject to nonoverlap.These steps are repeated for each particle in the simulation system before moving to a new pass.Whenever a trial move is successful, the effective volumes that are affected must be recomputed, making an impact to succeeding transitions [31].Periodic boundary conditions are assumed in the transitions to alleviate finite volume artifacts.A similar Monte Carlo procedure is employed in [23], which also incorporates between passes steps of particle insertion and removal, cluster translation, and parallel replica exchange to accelerate the convergence to equilibrium.Two examples of initial states for two different values of τ are given in Fig. 2. For the present Monte Carlo implementation of SHS, Fig. 3 indicates the accuracy of the PY approximation for the coordination number [Eq. ( 2)] in accordance with [25].

Formulation of the model
Overview.We use the commercial software PFC3D v5 by Itasca [26], which implements the original discrete element method described in [35].We simulate the confined uniaxial compression of a cohesive granular assembly of particles inside a cubic box of unit side length.The simulations are performed in the absence of gravity.Hence, a relatively high particle density ρ = 10 000 kg/m 3 (about that of silver) was chosen to improve the computational time since the time step is proportional to √ ρ.Initial conditions.The initial position and size of the particles are taken from the Monte Carlo procedure for different values of the volume fraction φ and stickiness τ .Note that the SHS configuration yields particles which are exactly at contact, which automatically leads to a mechanical equilibrium state in the absence of external loading for any new interaction force which is zero at contact.For each combination of (φ, τ ), simulations are conducted for three different realizations of initial configurations.The simulations were performed with N = 2048 particles, a number which was shown sufficient to avoid finite size effects.In addition, we verified that with N = 2048, the number of three realizations was sufficient as the system is already self-averaging (see Sec. IV D about finite size effects).Note that the radius of the particles is directly Boundary conditions.The granular samples were placed between two horizontal rigid walls of unit side length.The bottom wall is fixed while the top wall is translated at a fixed velocity v w = −0.01m/s.As will be shown, this velocity is small enough to consider the simulated mechanical behavior as quasistatic and independent of v w .Periodic boundary conditions were applied on the lateral walls.Hence, if a particle centroid falls outside of the model domain, it is translated back to the opposite side of the model.
Contact law.We used as cohesive contact law the PFC parallel bond model [36] also described in [37].The bonded part acts in parallel to the classical linear contact law [4,5].
For the unbonded, linear component [Fig.4(a), in gray], the normal force is the sum of a linear elastic and of a viscous term (spring-dashpot model), and the shear force is linear elastic with a Coulombian friction threshold.The mechanical parameters of this unbonded part of the model, namely, the contact elastic modulus E u and Poisson's ratio ν u , the restitution coefficient e u (viscous parameter), and the friction coefficient μ u are summarized in Table I.The value of the elastic modulus E u was chosen in such a way that the normal interpenetration at contacts are kept small, i.e., to work in the quasirigid grain limit [2,38] for the initial compression stage which we are interested in.However, this assumption does not stand during jamming (not studied here) due to extremely high normal forces.Concerning the normal restitution coefficient e u , we checked that the results presented below (macroscopic elastic modulus and strength) are uninfluenced by e u (in the range 0.1-0.9).This is due to the presence of the cohesive part of the contact law (see details below).Indeed, the restitution coefficient might have a strong influence for cases in which new contacts and collisions occur.In our case, the results are mostly driven by bond breaking, which explains why e has no influence on the results.However, the dense packing regime (jamming) which was not studied here might be influenced by the restitution coefficient.
The bonded part [Fig.4(a), in black] can be viewed as a cement joint with constant elastic modulus and Poisson's ratio E b and ν b acting at the contact.The bond failure behavior is characterized by its shear and tensile strength σ s and σ t .The maximum shear and tensile stresses σ max and τ max are calculated via beam theory [39] as follows: with N b and S b , the normal and shear forces in the bond, M b,1 the bending moment, M b,2 the twisting moment, r b the bond radius, A = πr 2 b the area of the bond, I = πr 4 b /4 the moment of inertia of the bond cross section, and J = πr 4  b /2 its the polar moment of inertia.The bond fails if σ t,max σ t or if σ s,max σ s (Fig. 4).If the bond fails in tension, the normal and shear forces are set to zero.If the bond fails in shear, the contact forces might not be altered but, only if the shear force does not exceed the friction limit and if the contact is in compression.The values of the mechanical properties used for the bonded part of the parallel bond model are given in Table I.Note that from the SHS equilibrium states, bonded contacts are created initially if the gap between particles is smaller than 1.5% of the particle diameter.Bonds were only introduced before the beginning of the simulation, and no new bonds were created during the simulations.Finally, as outlined in [40], our bond implementation, including bending, is similar to the effect of a rolling resistance (rolling friction coefficient) or particle shape and allows prevention of unphysical particle rotation.
Time stepping and elastic waves.The time step of the simulations was computed from particle properties according to t = √ m/k ≈ r √ ρ/E, where m, ρ, and r are the grain mass, density, and radius, and k and E the contact or bond stiffness and elastic modulus, respectively.The choice of this time step guarantees the algorithm stability [26,41].The magnitude of the elastic wave's celerity is c e ≈ √ E/ρ ≈ 30 m/s v w , ensuring the quasistatic assumption.

IV. MECHANICAL PROPERTIES OF COHESIVE SHS ASSEMBLIES A. Uniaxial compression
For an illustration of the generic behavior, Figs. 5 and 6 show examples for the evolution of the the compressive stress σ zz and cohesive coordination number z c with increasing axial deformation zz for different initial conditions (volume fraction φ and stickiness τ ).For low deformations, σ zz increases almost linearly with zz .During this initial phase, the number of cohesive contacts is constant and we refer to it as a quasielastic regime.This regime is characterized by a tangent modulus E, which is defined as the average slope in the σ zz − zz plane for which z c = z c,0 (no broken bonds).Due to the periodic boundary conditions, we expect this modulus to be in between the Young's modulus (obtained under unconfined compression) and the P-wave modulus (obtained under confined compression).For simplicity E is referred to as the elastic modulus henceforth.Note that the initial coordination number predicted by PY can be slightly different than that obtained with one realization of the SHS model, in particular, with low stickinesses (close to the coexistence line) as shown in Figs.4(b) and 5(b).After the quasielastic regime, a sudden drop of σ zz is observed corresponding to the failure.The failure corresponds to a small drop in the number of cohesive contacts, generally lower than 5% of the total number of cohesive contacts.The compressive stress at failure corresponding to the stress peak after the elastic phase is called the compressive strength σ c .After the failure and for zz 0.5, the stress gradually increases with deformation but is subject to fluctuations and sudden stress drops corresponding to consecutive bond failure events.For large deformations ( zz 0.5), σ zz increases strongly with zz in association with a strong decrease of the number of cohesive contacts.Overall, the compressive stress σ zz , elastic modulus E, and compressive strength σ c are increasing with increasing initial volume fraction φ 0 and decreasing with increasing stickiness τ .As a consequence [Fig. 1 and Eq. ( 2)], the latter quantities (σ zz , E, σ c ) are increasing with the initial cohesive coordination number z c and thus with the initial cohesive contact density ν c,0 = φ 0 z c,0 .
Although the macroscopic loading mode is compression, the introduction of cohesion and the loose character of the samples leads to large tensile and shear stresses in bonds, as shown in Fig. 7.The bonds are failing mostly due to tension or bending [Eq.( 4)] rather than shear or twisting [Eq.( 5)].Compressive stresses only appear for very high deformation ( zz 0.4) and thus high volume fractions.

B. Elastic modulus
We have analyzed the elastic modulus from the initial slope of the stress-strain curve and plotted it as a function of coordination number and contact density.The results are shown in Fig. 8.As predicted qualitatively in Sec.IV A, the elastic modulus increases with increasing cohesive coordination number z c,0 and with increasing volume fraction φ 0 .When plotted against the contact density ν c,0 = φ 0 z c,0 , all data points collapse on the same master curve that can be well described by a power law.By fitting the data, we obtain the regression  for the modulus of SHS-based cohesive assemblies.This leads to μ ≈ 4.9.To validate Eq. ( 6), we performed a detailed parametric study (in the Appendix) of the effect of the bond elastic modulus E b on the macroscopic modulus E, which confirms the linear relationship between E and E b .

C. Compressive failure
In the simulations, the compressive strength increases with both cohesive coordination number z c,0 and the initial volume fraction φ 0 [Fig.8(b)].Similar to the elastic modulus, when plotted against the cohesive contact density ν c,0 = φ 0 z c,0 , all σ c data points collapse on the same curve, which is best described by a power law.By fitting the data we obtain the regression for the compressive strength of SHS-based cohesive assemblies.Similar to the elastic modulus, we also confirmed (in the Appendix) that the macroscopic compressive strength σ c depended linearly on the bond tensile strength σ t .

D. Finite size effects
Before turning to the intermediate stage of compression, we provide some confidence that the scaling behavior found in Fig. 8 is not dominated by finite size effects.For a typical pair of values (φ = 0.3, τ = 0.13) (black square in Fig. 1) we have conducted simulations for different system sizes (total number N of spheres).We plotted the ratio E(ν c ,N )/E fit (ν c ) and σ c (ν c ,N )/σ c,fit (ν c ) for different N = 256, 512, 1024, 2048, 4096, 8092 in Fig. 9.Even though improved estimates would have been possible by choosing N = 4096, our choice of N = 2048 is reasonable and corresponds to a tradeoff between numerical effort and accuracy.

E. Plastic consolidation
After the initial elastic regime (σ zz < σ c ), the system turns to the intermediate stage of compression, which is generally characterized by the plastic consolidation curve describing irreversible compaction for an increasing stress [8,12].It is shown in Fig. 10 for different initial stickinesses and volume fractions.The plastic consolidation curve has the same characteristics as those reported in [12].For σ zz > σ c , irreversible plastic effects begin and we observe an almost logarithmic decrease of 1/φ.In this regime, volume fraction is related to the compressive stress according to where φ c = φ(σ c ) is the volume fraction of the system when the compressive stress σ zz meets the compressive strength σ c .The parameter λ is often called the plasticity index.We observe that the plasticity index is almost independent of the initial conditions (φ 0 , τ ) and that λ ≈ 0.5.The parameters φ(σ c ) and σ c , on the other hand, depend on the initial conditions as shown in Secs.IV B and IV C.

A. Deformation regimes in SHS-based cohesive materials
Similar to [12], the deformation of the SHS-based cohesive materials undergoes three characteristic stages with different mechanical behaviors.
First, for small deformations, no plastic (i.e., irreversible) effects occur and the samples behave like elastic solids.The elastic modulus scales with the initial contact density ν c,0 = φ 0 z c,0 according to a power law with exponent 4.9, independent of the bond size [cf.sensitivity analysis in the Appendix, Fig. 13(a)].This corroborates the results of [12], who highlighted a substantial increase of the elastic modulus with density and coordination number.This increase cannot be explained by classical theories [42] predicting a linear dependency of the elastic modulus with contact density.The previous models well reproduce the elastic behavior of cohesionless granular materials [12] for which the load is mostly carried by compressive and tangential contact forces.However, as shown in Fig. 7, cohesive systems mostly involve tensile and tangential force distributions in bonds, which is not accounted for in previously mentioned theories.Furthermore, note that in contrast to [12], the procedure to generate our samples based on the SHS prevents both inertial effects in the initial stage of compression and effects of the granular temperature on the initial coordination number, since the samples have initially no kinetic energy and the initial coordination number is prescribed [Eq.( 2)].
Second, after this elastic regime, a small number of bonds starts to break before the macroscopic plastic collapse of the samples.This collapse is characterized by a drop of the compressive stress σ zz (Figs. 5 and 6).The value of σ zz corresponding to this macroscopic failure is the compressive strength σ c of the sample.Similar to the elastic modulus, we showed that σ c scales with the initial contact density according to a power law with an exponent close to 3, also independent of the bond radius [cf.sensitivity analysis in the Appendix, Fig. 13(b)].As both the elastic modulus and compressive strength increase with increasing contact density ν c,0 , the amount of elastic deformation decreases with ν c,0 .After the macroscopic failure, the samples undergo a plastic consolidation regime [12].This regime is characterized by the plasticity index λ.In contrast with E and σ c , the plasticity index is independent of the initial conditions, i.e., the initial coordination number z c,0 and volume fraction φ 0 .This result is in line with the findings of [12], although the value of the plasticity index we find (λ ≈ 0.5) is higher than that of [12] (λ < 0.35).This discrepancy is certainly due to the difference of contact law used in the simulations.Our contact law does not allow for new cohesive contacts to be formed, in contrast to [12] for which new contacts can adhere.Hence, more plastic effects are expected in our case and thus a higher value of λ.During this phase, the bond normal and shear stress distributions [Eqs.( 4) and ( 5)] widen as axial deformation increases (Fig. 7).The normal stress is essentially in tension, which is the main mode of failure of cohesive bonds.Indeed, shear stresses in the bonds are approximately 50% lower than tensile ones.This result confirms the findings of [12], which argued that the deformation of loose systems was mostly controlled by the bending of thin arms which connect denser zones moving like rigid bodies.The plastic consolidation regime is commonly discussed as the plateau region in collapsible, cellular materials where strain localization and the propagation of crushing bands are predicted [43] if the degree of structural heterogeneity is sufficiently small.We did not observe these features for our low-density SHS assemblies; a systematic analysis is, however, beyond the scope of the present work.
Finally, after the latter plastic consolidation regime, the onset of jamming occurs and the compressive stress increases strongly while the volume fraction and axial deformation slightly increase.The granular material becomes extremely rigid and the mechanical behavior is independent of the initial conditions.In this regime, the few remaining bonds have stress states not only in tension but also in compression, as shown in Fig. 7 for an axial deformation zz = 0.7.Note that we did not intend to analyze the details of the jamming transitions as we are mainly interested in the two first regimes of deformation.In particular, for this very high pressure regime, the quasirigid assumption of our particles is no longer relevant.In addition, the monodisperse size distribution associated with the SHS generation method induces crystallization effects.
Our analysis of loose and cohesive granular materials under uniaxial compression represents a first step towards a more complete characterization under mixed stress states (shear and compression).This will allow in the future to evidence the influence of contact density on the yield surface and make a connection to critical state soil mechanics [8,44,45].

B. The role of contact density and scaling exponents
Granular materials are commonly analyzed in terms of volume fraction and coordination number as key microstructural descriptors.In quasirigid, bonded sphere assemblies the scaling of E and σ c cannot be determined by volume fraction alone since the elastic response is predominantly caused by the cohesive bonds and their spatial distribution.This is reflected by classical considerations from the engineering literature, where under the assumption of an affine deformation field, macroscopic strains can be translated directly to the strains in the cohesive bonds.This implies that the elastic modulus E, which has dimensions of strain energy per unit volume, should scale linearly with the the number of bonds per unit volume, i.e., with the contact density ν c = z c φ via E ≈ ν μ c (9) with μ = 1, which is exactly the classical result of Walton [42].The relevance of contact density ν c was recently confirmed by [18] and [7], which showed a significant influence on the elastic moduli.In particular, [18] suggested that ν c should be used as indicator of the internal state of granular packings.However our initial states bear spatial fluctuations in both particle volume fraction and the coordination number of cohesive bonds, leading to strong deviations from the prediction μ = 1 for a perfectly homogeneous configuration.The presence of heterogeneities is also reflected by the stress distribution (Fig. 7), which is far from being sharply peaked for all SHS parameters.Spatial fluctuations of the mechanically active units (cohesive bonds) induce heterogeneities in the strain field [46], and deviations from Eq. ( 9) with μ = 1 should be expected.
Similar arguments apply when the system fails macroscopically.Under the assumption of an affine deformation field, bond failure directly dictates the macroscopic behavior.This is captured by the Rumpf model [47], which can be recast into a scaling with contact density according to σ c ≈ ν α c (10) with α = 1, as the analog of Eq. ( 9) for the compressive strength, again oversimplifying the observed behavior.
To discuss the value of the observed exponents μ,α it is natural to consider related nongranular systems.The scaling (Fig. 8) is expected to originate from the underlying structural transitions, namely, percolation and phase separation, which dominate the structure of SHS similar to short-range attractive colloids [23,48].The existence of a percolation transition [22] in the Baxter model would suggest a scaling E ≈ (φ − φ crit ) μ in the vicinity of the percolation line τ perc in Fig. 1.The careful analysis of the SHS phase diagram [25], however, shows the poor performance of PY for the percolation and coexistence line.Thus the precise position of our initial states (chosen according to PY, Fig. 1) in the true phase diagram is difficult to assess.Combining both dependencies of the moduli shown in Fig. 8, our results would imply a power law E ≈ φ 4.9 0 and σ c ≈ φ 3 0 if regarded only as a function of initial volume fraction φ 0 .It is often observed [49] that gelling systems with short-range attractions exhibit power-law behavior E ≈ φ μ of the elastic moduli, even without reference to the onset of rigidity φ crit .The elastic exponent commonly varies between 3 < μ < 5 [49], consistent with our finding μ ≈ 4.9.Numerous scaling arguments have been proposed for the elasticity exponent in gelling systems in the previous decades.The value of μ thereby reflects the connectivity of the percolating cluster and details of the pair interactions (or contact law) [50].The cubic dependence of the elastic modulus of individual samples on the bond radius r b (cf.sensitivity analysis in the Appendix, Fig. 11) reflects the primary influence of the bond-bending term in the interaction law [Eq.( 4)].The observed exponents are independent of r b , at least in the range of bond radii considered in the sensitivity analysis (Fig. 13).Recent work [51] shows that it is experimentally feasible to fine tune bond parameters via "crosslinking" of glass beads.For these systems, elasticity and failure under uniaxial compression is mainly controlled by bond properties, even for higher volume fractions.We expect that the scaling of the mechanical properties with contact density observed in the present study is likely a robust feature and may be found for other assembling methods as well.
The value of the strength exponent α is typically smaller than μ [49], also consistent with our findings.The failure point can be inferred from the limit of linearity for the strain (see Ref. [49], and references therein) which separates reversible from irreversible deformations.We checked that our moduli E and σ c are estimated in a stage where compressive strains are reversible (up to a few percent).

VI. CONCLUSIONS
Three-dimensional uniaxial compression simulations of loose and cohesive granular assemblies were performed using the discrete element method (DEM).Initial configurations were obtained from Baxter's sticky hard spheres (SHS) model, which allows one to achieve mechanically rigid configurations at very low volume fractions.A unique feature of the SHS model is to generate initial states with different microstructures for identical volume fractions φ but with different coordination numbers z c , offering two independent controls on the contact density z c φ as the key quantity.
In the DEM simulations, samples undergo three deformation regimes: (i) the quasielastic regime for low deformations and pressures, for which almost no bond-breaking events occur; (ii) the plastic consolidation regime; and (iii) the dense packing regime corresponding to the jamming transition.We mainly focused on the analysis of the initial elastic modulus E and compressive failure strength σ c .Our results highlight the universal control of the cohesive contact density ν c = z c φ on the mechanics in regime (i), where E and σ c follow a power law as a function of ν c .On the other hand, the plasticity index λ, characterizing the plastic consolidation curve after failure in regime (ii), is independent of the initial conditions.
We believe that other assembling procedures would come to a similar scaling picture with different exponents which reflect the connectivity of the initial, mechanically active backbone.One advantage of using SHS is its wide use, which allows one to compare the mechanics of cohesive granular assemblies with other heterogeneous materials.Another advantage is available analytical approximations, e.g., for the pair correlation function which allows mapping of (bi-)continuous two-phase materials onto discrete SHS sphere assemblies by matching correlation functions from x-ray tomography via parameter optimization (radius, stickiness, or contact density) [52].Accordingly, our procedure of SHS-based generation of DEM initial states opens new perspectives for the understanding and mechanical modeling of other loose and cohesive materials, such as snow, rocks, soils, cement, asphalts, etc., with potential applications in natural hazards and engineering.simulations presented in the main text were performed for a ratio of tensile to shear strength σ t /σ s = 1.The effect of σ t /σ s is shown in Fig. 12(b).The macroscopic compressive strength σ c decreases with increasing σ t /σ s .For σ t /σ s 2, the failure occurs mostly because of tensile failure in bonds, while for σ t /σ s 2, the failure occurs due to shear and tensile failure of the bonds.In line with the effect of the bond radius r b on E, the macroscopic compressive strength σ c increases with ∼r 3 b .

Scaling laws
Finally, the effect of bond-bending forces on the two scaling laws obtained for the elastic modulus E [Eq.(7)] and compressive strength σ c [Eq. ( 9)] was studied.The bond radius was thus modified for four simulation points with different initial volume fraction φ 0 and stickinesses τ .Figure 13 shows the result of this analysis.Each point represents the average value of three different realizations of the initial state.The analysis reveals that the existence of a power-law relationship and the power-law exponent does not depend on r b /r p and is thus independent of the intensity of bond-bending forces, which strongly depend on r b [Eqs.( 4) and ( 5 (7)] and (b) the normalized compressive strength σ c /σ t [Eq.(9)].The additional simulations were performed for (φ 0 , τ )={(0.2,0.14) (0.3, 0.13) (0.35, 0.07) (0.35, 0.52)}, except for the effect of contact density on the compressive strength for r b /r p = 0.25 for which the simulations were performed for (φ 0 , τ )={(0.3,0.11) (0.3, 0.13) (0.35, 0.07) (0.35, 0.14)}.

FIG. 1 .
FIG. 1. Schematic of the SHS phase diagram in the φ,τ plane together with PY results and the choice of the simulation points.

FIG. 3 .
FIG. 3. Coordination number z SHS c evaluated from the Monte Carlo configurations vs the PY approximation z PY c .

FIG. 4 .
FIG. 4. Representation of the PFC parallel bond model used in the simulations.The bonded part is represented in black while the unbonded part is represented in gray.(b) Bond normal force N b as a function of the normal interpenetration δ n scaled by the bond radius r b .(c) Bond shear force ||S b || as a function of tangential interpenetration δ s scaled by the bond radius r b .(d) Bond-bending moment ||M b,1 || as a function of bending rotation θ 1 scaled by the bond radius r b .(e) Twisting moment ||M b,2 || as a function of twist rotation θ 2 scaled by the bond radius r b .

21 FIG. 5 .
FIG. 5. (a) Compressive stress σ zz vs axial deformation zz for φ 0 = 0.3 and for one realization of the initial state.(b) Cohesive coordination number z c vs axial deformation for φ 0 = 0.3 and for one realization of the initial state.The dotted lines correspond to the PY approximation for the initial state.

35 FIG. 6 .
FIG. 6.(a) Compressive stress σ zz vs axial deformation zz for τ ∈ [0.1 0.11] and for one realization of the initial state.(b) Cohesive coordination number z c vs axial deformation for τ ∈ [0.1 0.11] and for one realization of the initial state.The dotted lines correspond to the PY approximation for the initial state.

7 FIG. 7 .
FIG.7.Probability distribution function of the normal stress in bonded contacts σ t,max (positive in tension) normalized by the bond tensile strength σ t vs σ t,max /σ t .Inset: Probability distribution function of the shear stress in bonded contacts σ s,max normalized by the bond shear strength σ s vs σ s,max /σ s .The simulation was performed for φ 0 = 0.3 and τ = 0.11 for different deformation stages.

FIG. 12 .
FIG. 12. Dependence of the compressive strength σ c with (a) the bond tensile strength σ t ; (b) the ratio between tensile and shear strengths of the bond σ t /σ s ; and (c) the bond relative radius r b /r p .The simulations were performed for φ 0 = 0.3 and τ = 0.13.

TABLE I .
Mechanical parameters used in the simulations for the parallel bond model: E u contact elastic modulus, ν u contact Poisson's ratio, μ u contact intergranular friction, e u contact normal restitution coefficient, E b bond elastic modulus, ν b the bond Poisson's ratio, σ s bond shear strength, σ t bond tensile strength, r b bond radius, and r particle radius.
determined by the volume fraction φ and number of particles N in a box of unit size.