Microstructure evolution of compressible granular systems under large deformations
We report three-dimensional particle mechanics static calculations that predict the microstructure evolution during die-compaction of elastic spherical particles up to relative densities close to one. We employ a nonlocal contact formulation that remains predictive at high levels of confinement by removing the classical assumption that contacts between particles are formulated locally as independent pair-interactions. The approach demonstrates that the coordination number depends on the level of compressibility, i.e., on the Poisson’s ratio, of the particles. Results also reveal that distributions of contact forces between particles and between particles and walls, although similar at jamming onset, are very different at full compaction. Particle-wall forces are in remarkable agreement with experimental measurements reported in the literature, providing a unifying framework for bridging experimental boundary observations with bulk behavior.
] ] ]
The microstructure of confined granular media is typically inhomogeneous, anisotropic and disordered. Under external loading, these systems exhibit a non-equilibrium jamming transition from a liquid-like to a solid-like state [28, 24, 49, 39, 45]. Under increasing confinement, these amorphous solids support stress by spatial rearrangement and deformation of particles and by the development of inhomogeneous force networks [30, 17, 32, 3, 29, 41, 11]. Understanding and predicting the formation and evolution of such microstructure under finite macroscopic deformations and at high levels of confinement remains a fundamental goal of granular mechanics [18, 10, 42, 4, 13, 40].
In this paper, we report three-dimensional particle mechanics static calculations that enable us to predict microstructure evolution during die-compaction of elastic spherical particles up to relative densities close to one. We employ a nonlocal contact formulation  that remains predictive at high levels of confinement by removing the classical assumption that contacts between particles are formulated locally as independent pair-interactions. We specifically study a noncohesive frictionless granular system comprised by 40,000 weightless elastic spherical particles with radius mm, Young’s modulus GPa, and Poisson’s ratio . The granular bed, which is numerically generated by means of a ballistic deposition technique , is constrained by a rigid cylindrical container of diameter mm. The systems is deformed under die-compaction up to a relative density of . Assuming a sufficiently small compaction speed, we consider rate-independent material behavior and we neglect traveling waves, or any other dynamic effect . The deformation process is therefore described by a sequence of static equilibrium configurations. In this work we employ 125 quasi-static load steps and we consider two materials that only differ in their level of compressibility, namely and .
The paper is organized as follows. The particle mechanics approach used to generate a sequence of static equilibrium configurations of granular systems at high levels of confinement is presented in Section 2. The evolution of microstructural statistical features of these equilibrium configurations is investigated in Section 3. Specifically, we study the evolution of the mechanical coordination number (number of non-zero contact forces between a particle and its neighbors) in section 3.1, punch and die-wall pressures in section 3.2, the network of contact forces in section 3.3, and the network of contact radiuses in section 3.4. Finally, a summary and concluding remarks are collected in Section 4.
2 Particle mechanics approach to granular systems at high levels of confinement
Due to the high level of confinement experienced by the particles in the system, we adopt a nonlocal contact formulation  that accounts for the interplay of deformations due to multiple contact forces acting on each single particle. This nonlocal formulation removes the classical assumption that contacts between particles are formulated locally as independent pair-interactions. Its predictions are in remarkable agreement with detailed finite-element simulations of a neo-Hookean solid extended to the compressible range and with experimental observations of a rubber sphere that exhibits no permanent deformations. This good agreement at moderate to large levels of confinement, however, is not obtained with Hertz theory (see Figure 1 and  for further details). It bears emphasis that the nonlocal contact formulation is not an extension of Hertz theory to finite deformations but rather a systematic approach for accounting the interaction between contact interfaces. In a manner analogous to interacting cracks or dislocations in elastic media, the nonlinear behavior is confined to small regions and the interaction between these regions is through an elastic field of deformations.
An equilibrium configuration is defined by the solution of two sets of coupled nonlinear equations. The first set of equations corresponds to static equilibrium of the granular system, that is sum of all elastic contact forces acting on each particle equals zero. The second set of equations accounts for the contribution to each contact interface of the nonlocal mesoscopic deformations induced by all other contact forces acting on the particles, that is for the nonlocal terms . Thus, the nonlinear system of equations is given by
where and are the position and all the neighbors of particle , respectively, by definition, and . In the above system of equations, and correspond to Hertz theory and to the nonlocal contact formulation , respectively, with
where is the angle between the coordinate points of particles , and (see Figure 2).
A sequential strategy is proposed to treat the coupled problem (see Algorithm 1). The equations of static equilibrium are solved for , assuming fixed , by employing a trust-region method [15, 16] that successfully overcomes the characteristic ill-posedness of the problem (e.g., due to metastability ). The basic trust-region algorithm requires the solution of a minimization problem to determine the step between iterations, namely the trust-region step. This minimization problem is of the form , where is the trust-region step, is a trust-region radius and is a quadratic function that represents a local model of the objective function about , that is
with and the global force vector and stiffness matrix at —the first term in the above equation is not required in the minimization problem. It is worth noting that the trust-region step is not necessarily in the direction of a quasi-Newton step and that the trust-region radius acts as a regularization term that controls the growth in the size of the least squares solution observed in most ill-posed . Trading accuracy for performance, Byrd , among others, proposed to approximate the minimization problem by restricting the problem to a two-dimensional subspace. Furthermore, the two-dimensional subspace can be determined by a preconditioned conjugate gradient process and the trust-region radius can be adjusted over the iterative process (see, e.g., ). Here we adopt the implementation available in MATLAB R2014a Optimization Toolbox.
The nonlocal terms are solved by a few fixed point iterations of the second set of equations in (2), assuming a given (see  and Algorithm 2). It is worth noting that an equilibrium configuration is not obtained by artificially dumped or cooled-down dynamic processes but rather by iterative solvers that follow the energy landscape around the solution of static equilibrium. Each iteration of the trust-region method involves the solution of (subject to a given tolerance TOL2) and the evaluation of an approximate solution for . A notable feature of the nonlocal contact formulation is that, when nonlocal effects are neglected, it reduces to Hertz theory. Under such simplification, an equilibrium configuration is readily obtained from the solution of the first set of equations in (2), assuming .
3 Post-jamming behavior of compressible granular systems
We next study the post-jamming behavior of two granular systems that only differ in the level of compressibility of the particles (i.e., in ). The geometric configuration of the initial granular bed is the same for both systems, and simulations are performed for both Hertz theory and the nonlocal contact formulation (each simulation takes several months of CPU time). Deformations at the particle scale are within the validated range of the nonlocal contact formulation, that is smaller than 40 (see Figure 1), for relative densities smaller than 0.90. We investigate the evolution of statistical features of the mechanical coordination number (number of non-zero contact forces between a particle and its neighbors), punch and die-wall pressures, the network of contact forces, and the network of contact radiuses. Error bounds in tables and plots represent the 95% confidence envelope for fitted coefficients.
3.1 Mechanical coordination number
The mean coordination number evolves as a power law of the following form
where is the critical relative density, is the minimal average coordination number and is the critical exponent. This well-known critical-like behavior has an exponent consistent with for different pair-interaction contact laws, polydispercity and dimensionality of the problem [38, 39, 17]. Figure 3(a) shows the results obtained from the particle contact mechanics simulations for Hertz theory and their best fit to equation (3). For , jamming occurs at and with . The fit to numerical results is excellent not only near jamming but also at large relative densities [see Figure 3(b)]. It is worth noting that the isostatic condition for frictionless packings implies a critical coordination number equal to 6 and a critical density close to 0.64. Our initial powder bed, however, is generated such that a subset of particles is jammed, namely the backbone, while the remainder of the particles are not jammed but locally imprisoned by their neighbors, referred to as rattlers. In this work, we do not exclude rattlers from packings, when calculating packing densities, mean coordination numbers and networks of contact forces. Therefore, compaction is initially characterized by rearrangement and very small deformations of the particles, and it is followed by an early jamming transition at , reaching the isostatic coordination number at a relative density very close to 0.60. There exists a body of work that indicates that depends on the protocol used for obtaining jammed configurations and that monodisperse systems are prone to crystallization [2, 12, 43, 51]. While a detailed analysis concerning these points is a worthwhile direction of future work, specially the investigation of polydisperse packings for which some results are observers to be independent of the generation protocol, here we restrict our discussion to post-jamming behavior and to one preparation protocol. Results for are only marginally different.
Hertz theory, and any other contact law formulated locally as pair-interactions, does not account for the interplay of deformations due to multiple contact forces acting on a single particle. As relative density increases, these nonlocal deformations are responsible for the formation and breakage of contacts between particles which in turn manifest as a crossover in the local coordination number [see top insert in Figure 3(a) for two bi-dimensional examples]. This local crossover is predicted by the nonlocal contact formulation and it emerges in the evolution of the mean mechanical coordination number depicted in Figure 3(a). For highly compressible materials (), the evolution of is below the predictions of Hertz theory, i.e., some contacts are broken as the level of confinement increases. In sharp contrast, for nearly incompressible materials (), the evolution of is above the predictions of Hertz theory, i.e., new contacts are formed due to volume conservation at the particle scale. Therefore, while equation (3) is valid and independent of any material property near jamming, its validity fails at high levels of confinement where strongly depends on .
The probability distribution of the local coordination number is best described by the granocentric model proposed in [14, 37]. Here we restrict our attention to monodisperse systems and approximate by a binomial distributed variable characterized by its mean and variance, that is
where represents the mean number of (closest) neighbors and corresponds to the probability of finding a contact among these neighbors. Inserts in Figure 3(a) show that this is an accurate approximation for both formulations at small and large relative densities—though distributions for the nonlocal contact formulation are not shown in the figure. The numerical results also suggest the existence of an universal relationship between and of the following form
with and , for both formulations and Poisson’s ratios [see Figure 3(b) for exact numerical values obtained for hertzian interactions].
Finally, we investigate the sensitivity of these results to the numerical tolerance used to identify the non-zero contact forces that contribute to the mechanical coordination number. We thus define as the inter-particle deformation above which a contact force is identified as non-zero. We also investigate how the ordering effect of the container walls penetrates into the bulk of the powder bed. We thus define as the distance between the boundary and the bulk region. Figure 4(a) shows that the evolution of the mean coordination number is insensitive to the choice of at high relative densities and barely sensitive at small densities, for values of equal to , , and and for . Figure 4(b) indicates that the effect of the boundary on the mean coordination number only penetrates one or two particle radiuses into the bulk, for equal to . Therefore, the results presented in Figure 3 correspond to and .
3.2 Applied pressure
The pressures applied by the punches and the reaction at the die wall are effective macroscopic variables predicted by the particle contact mechanics simulation. These predictions are presented in Figure 5(a) and follow a power law of the following form
where the pressure is evaluated at the deformed configuration, is obtained from the evolution of , the coefficients in the first term of the equation (i.e., and ) are best-fitted to the numerical results for Hertz theory, and the remaining coefficients (i.e., and ) are best-fitted to predictions of the nonlocal contact formulation [see Figure 5(b) for numerical values]. It is interesting to note that stiffnesses and critical exponents for hertzian systems are well approximated by
near jamming (see, e.g., [38, 17, 32, 33]) and at large relative densities. Also, and are consistent with the exponents in a series expansion of the pressure applied in the configuration presented in Figure 1 (see equation 12 in ). An analytical derivation of is a worthwhile direction of future work, for which one may expect to only depend on and thus to be in agreement with studies presented in  for granular crystals or highly packed granular lattices. In a similar manner, one may expect that the behavior observed in Figure 3 for solely depends on . A detailed analysis concerning this point is however beyond the scope of this work.
3.3 Network of contact forces
The probability distribution of the contact forces is assumed to have an algebraic tail with exponent and to be translated by , that is
where is the non-dimensional contact force and is a shape parameter that controls the behavior of weak forces—for , is finite if ; and for , is finite. Contact forces between particles and between particles and walls are studied separately and non-dimensionalized by their own average value . Figure 6 shows that, as relative density increases, the tails of both probability distributions transition from exponential () to Gaussian (). However, the distributions, although similar at jamming onset, are very different at full compaction. The particle-wall force distribution becomes increasingly symmetric, that is its mode approaches its mean value. In sharp contrast, the particle-particle force distribution shows a peak at small forces and a local maximum, or a plateau, whose location saturates to . These results agree with previous studies that noted a transition from exponential to Gaussian tails [1, 6, 7, 19, 25, 32, 47, 50]. They also extend to three-dimensional systems the observation that distributions measured at the walls do not reflect accurately distributions of the bulk in general (see, e.g., [44, 48] for two-dimensional studies).
A best fit of the numerical results suggests the following further simplification
Thus, after best fitting the numerical results of the particle contact mechanics simulations to (8) with (9-10), it is evident from Figure 6 that two parameters fully-describe each probability distribution with remarkable accuracy. Furthermore, the numerical results are in excellent agreement with experimental measurements reported in the literature for one-dimensional compaction of spherical deformable particles in cylindrical containers. These measurements are obtained at, or near, the wall using a carbon paper technique [19, 36, 32], or confocal microscopy , respectively—with the remarkable exception of , and recently , which are obtained within the bulk of monodisperse emulsions by confocal microscopy. Naturally, the statistical significance obtained for particle-wall interactions is smaller than the one obtained for particle-particle interactions—e.g., in this study, the ratio between the number of each type of interaction remains almost constant during compaction and equal to 3%.
Figure 7 shows the evolution of for particle-particle interactions and of for particle-wall interactions. Interestingly, the figure suggests a critical-like behavior for the parameter with critical exponent and critical density . It also reveals that predictions of the nonlocal contact formulation depart from those obtained with Hertz theory at moderate to high levels of confinements (e.g., saturates at 1.8, not reaching truly Gaussian values). The systematic investigation of this critical-like behavior, and the elucidation of the material properties and microstructural features that control it, are worthwhile directions of future research. In addition, the figure revelas that is larger than 0.6 and saturates at 0.80, whereas is found to be 1. This is in agreement with recent theoretical and numerical calculations (see, e.g., [5, 27] and references therein). It is worth noting that the probability distribution of particle-wall contact forces thus takes a finite value at zero proportional to . Finite values for have been reported in  for particle-particle forces.
Finally, Figure 8 illustrates the network of forces in the three-dimensional granular system for four different relative densities. We only visualize non-dimensional contact forces larger than 2.5, that is the forces that carry most of the load in the system. It is clearly observed in the figure that at low relative density or confinement the network is inhomogeneous in space and characterize by a broad range of forces. In contrast, at high relative density or confinement the network becomes more homogeneous in space and the range of observed forces is narrower . These observations are consisten with the transition from exponential to Gaussian tails in the distribution of particle-particle interactions under increasing relative density.
3.4 Network of contact radiuses and areas
The probability distribution of the contact radiuses is assumed to have an algebraic tail with exponent and to be translated by , that is
where we have changed variables in (8) using . We note that the hertzian relationship does not hold for the average values (i.e., ) and thus the non-dimensional values are related by the following inequality . The identification of tighter bounds for the inequality and how these bounds can be used to bound are however beyond the scope of this work. Figure 9 shows that, as relative density increases, the range of observed contact radiuses is narrower. This transition however occurs faster than in the case of contact forces for particle-particle interactions (cf. Figure 6). In accordance with the behavior of the network of forces, the figure indicates that the distributions, though similar at jamming onset, are very different at full compaction. Figure 9 also shows a best fit of the numerical results to equation (11) for reference purposes. However, the simplifications proposed for the network of forces, i.e., equations (9-10), do not hold true for the network of contact radiuses. It is work noting that the choice of does influence the tail of the probability distribution at when the system is at jamming onset, see Figure 9(a).
4 Summary and discussion
We have reported three-dimensional particle mechanics static calculations that predict the microstructure evolution of monodisperse elastic spherical particles during die-compaction up to relative densities close to one. We employ a nonlocal contact formulation that remains predictive at high levels of confinement by removing the classical assumption that contacts between particles are formulated locally as independent pair-interactions. Our approach demonstrates that the coordination number depends on the level of compressibility, i.e., on the Poisson’s ratio, of the particles and thus its scaling behavior is not independent of material properties as previously thought. Our results also reveal that distributions of contact forces between particles and between particles and walls, although similar at jamming onset, are very different at full compaction. Both distribution tails transition from exponential to Gaussian under increasing relative density—with particle-wall forces exhibiting an earlier transition. Particle-wall forces are in remarkable agreement with experimental measurements reported in the literature, providing a unifying framework for bridging experimental boundary observations with bulk behavior. The evolution of these microstructural features is not otherwise amenable to direct experimental determination. Furthermore, prediction of such evolution is key to better design, optimize and control many manufacturing processes widely used in pharmaceutical, energy, food, ceramic and metallurgical industries.
We close by pointing out some limitations of our approach and possible avenues for extensions of the analysis.
First, we have generated the initial powder bed by means of a ballistic deposition technique and progressively moving the upper punch until particle rearrangement and small deformations lead to a jammed configuration. This is not the only—perhaps even the best—available protocol for obtaining jammed configurations in rigid containers. A case in point consists in modeling the filling process (see, e.g., [31, 53]) by solving the equations of motion of the Lagrangian system of particles with a numerical time integrator (see, e.g., [20, 22, 23] and references there in). However, this approach, though more realistic, may also result in a jammed backbone and a number of rattlers. The systematic investigation of the effect of the filling protocol on the formation and evolution of the microstructural features studied here at high relative densities is a worthwhile direction of future research.
Second, we have restricted attention to die-compaction of monodisperse system constrained by a rigid cylindrical container and two flat punches.The systematic extension of this analysis to dies and punches of other shapes, hydrostatic compaction conditions, polydispersity and mixtures of particles with different material properties are interesting topics to study.
Third, although our analysis has focused on coordination number, punch and die-wall pressures, and the network of contact forces and radiuses, there are other microstructural features of theoretical and practical importance. For example, it is relevant to study the evolution of material-fabric tensor, force-fabric tensor and pore size distribution during die compaction.
The authors gratefully acknowledge the support received from the NSF ERC grant number EEC-0540855, ERC for Structured Organic Particulate Systems. MG also acknowledges support from Purdue University’s startup funds.
-  Antony S.J. Evolution of force distribution in three-dimensional granular media, Physical Review E 2000; 63:1, 011302.
-  Baranau V. and Tallarek U. On the jamming phase diagram for frictionless hard-sphere packings, Soft matter 2014; 10, 7838.
-  Blair D.L., Mueggenburg N.W., Marshall A.H., Jaeger H.M., and Nagel S.R. Force distributions in three-dimensional granular assemblies: Effects of packing order and interparticle friction, Physical Review E 2001; 63:4, 041304.
-  Blumenfeld R. and Edwards S.F. On granular stress statistics: Compactivity, angoricity, and some open issues, Journal of Physical Chemistry B 2009; 113, 3981–3987.
-  Bo L., Mari R., Song M., and Makse H.A. Cavity method for force transmission in jammed disordered packings of hard particles, Soft Matter 2014; 10, 7379–7392.
-  Brujić, J., Edwards S.F., Grinev D.V., Hopkinson I., Brujić D., and Makse H.A. 3D bulk measurements of the force distribution in a compressed emulsion system, Faraday Discuss 2003; 123, 207–220.
-  Brukić J., Edwards S.F., Hopkinson I., and Makse H.A. Measuring the distribution of interdroplet forces in a compressedemulsion system, Physica A 2003; 327, 201–212.
-  Brodu N., Dijksman J.A., and Behringer R.P. Spanning the scales of granular materials: Microscopic force imaging, Nature Communications 2015; 6, 6361.
-  Byrd R.H., Schnabel R.B., and Shultz G.A. Approximate solution of the trust region problem by minimization over two-dimensional subspaces, Mathematical Programming 1988; 40, 247–263.
-  Cates M.E., Wittmer J.P., Bouchaud, J.-P., and Claudin P. Jamming, force chains, and fragile matter, Physical Review Letters 1998; 81:9, 1841–1844.
-  Chan S.H. and Ngan A.H.W. Statistical distribution of contact forces in packings of deformable spheres, Mechanics of Materials 2005; 37, 493–506.
-  Chaudhuri P., Berthier L., and Sastry S. Jamming transitions in amorphous packings of frictionless spheres occur over a continuous range of volume fractions, Physical Review Letters 2010; 104, 165701.
-  Ciamarra M.P., Richard P., Schröter M., and Tighe B.P. Statistical mechanics for static granular media: open questions, Soft Matter 2012; 8, 9731–9737.
-  Clusel M., Corwin E.I., Siemens A.O.N, and Brujic J. A granocentric model for random packing of jammed emulsions, Nature 2009; 460, 611–615.
-  Coleman T.F. and Li Y. An interior trust region approach for nonlinear minimization subject to bounds, SIAM Journal on Optimization 1996; 6:2, 418?-445.
-  Conn, A. and Gould, N. and Toint, P. Trust Region Methods, Society for Industrial and Applied Mathematics, 2000.
-  Durian D.J. Foam mechanics at the bubble scale, Physical Review Letters 1995; 75:26, 4780 -4783.
-  Edwards S.F. and Oakeshott R.B.S. Theory of powders, Physica A 1989; 157, 1080–1090.
-  Erikson J.M., Mueggenburg N.W., Jaeger H.M., and Nagel S.R. Force distributions in three-dimensional compressible granular packs, Physical Review E 2002; 66, 040301(R).
-  Gonzalez M., Schmidt B. and Ortiz M. Energy-stepping integrators in Lagrangian mechanics, International. Journal for Numerical Methods in Engineering 2010; 82:2, 205–241.
-  Gonzalez M., Yang J., Daraio C., and Ortiz M. Mesoscopic approach to granular crystal dynamics Physical Review E 2012, 85, 016604.
-  Gonzalez M. and Cuitiño A.M. A nonlocal contact formulation for confined granular systems, Journal of the Mechanics and Physics of Solids 2012; 60, 333–350.
-  Hairer E, Lubich C, Wanner G. Geometric numerical integration. Structure-Preserving Algorithms for Ordinary Differential Equations Springer Series in Computational Mathematics (2nd edn), vol. 31. Springer: Berlin, 2006.
-  Jin Y. and Makse H.A. A first-order phase transition defines the random close packing of hard spheres, Physica A 2010; 389:23, 5362 –5379.
-  Jorjadze I., Pontani L.-L., and Brujic J. Microscopic Approach to the Nonlinear Elasticity of Compressed Emulsions, Physical Review Letters 2013; 110:4, 048302.
-  Jullien R. and Meakin P. Simple-models for the restructuring of 3-dimensional ballistic aggregates. Journal of Colloid and Interface Science 1989; 127:1, 265– 272.
-  Lerner E., During G., and Wyart M. Low-energy non-linear excitations in sphere packings, Soft Matter 2013; 9, 8252–8263.
-  Liu A.J. and Nagel S.R. Nonlinear dynamics: Jamming is not just cool any more, Nature 1998; 396:6706, 21–22.
-  Liu C.-H., Nagel S.R., Schecter D.A., Coppersmith S.N., Majumdar S., Narayan O., and Witten T.A.. Force fluctuations in bead packs, Science 1995; 269, 513–513.
-  Majmudar T.S. and Behringer R.P. Contact force measurements and stress-induced anisotropy in granular materials, Nature 2005; 435, 1079–1082.
-  Mateo-Ortiz D., Muzzio F.J. and Méndez R. Particle size segregation promoted by powder flow in confined space: The die filling process case, Powder Technology 2014; 262, 215–222.
-  Makse H.A., Johnson D.L., and Schwartz L.M. Packing of compressible granular materials, Physical Review Letters 2000; 84:18, 4160–4163.
-  Mason T.G., Lacasse M.-D., Grest G.S., Levine D., Bibette J., and Weitz D.A. Osmotic pressure and viscoelastic shear moduli of concentrated emulsions, Physical Review E 1997; 56:3, 3150 -3166.
-  Mehta A. Granular physics, Cambridge University Press, Cambridge, 2007.
-  Moré J.J. and Sorensen D.C. Computing a trust region step, SIAM Journal on Scientific and Statistical Computing 1983; 4:3, 553?-572.
-  Mueth D.M., Jaeger H.M., and Nagel S.R. Force distribution in a granular medium, Physical Review E 1998; 57:3, 3164–3169.
-  Newhall K.A., Jorjadze I., Vanden-Eijnden E., and Brujic J. A statistical mechanics framework captures the packing of monodisperse particles Soft Matter, 2011; 7, 11518.
-  O Hern C.S., Langer S.A., Liu A.J., and Nagel S.R. Random packings of frictionless particles, Physical Review Letters 2002; 88:7, 075507.
-  O Hern C.S., Silbert L.E., Liu A.J., and Nagel S.R. Jamming at zero temperature and zero applied stress: The epitome of disorder, Physical Review E 2003; 68:1, 011306.
-  Ostojic S., Somfai E., and Nienhuis B. Scale invariance and universality of force networks in static granular matter Nature 2006; 439, 828–830.
-  Radjai F., Jean M., Moreau J.-J., and Roux S. Force Distributions in Dense Two-Dimensional Granular Systems, Physical Review Letters 1996; 77:2, 274–277.
-  Richard P., Nicodemi M., Delannay R., Ribiere P., and Bideau D. Slow relaxation and compaction of granular systems. Nature Materials 2005; 4:2, 121–128.
-  Schreck C.F., O?Hern C.S., and Silbert L.E. Tuning jammed frictionless disk packings from isostatic to hyperstatic Physical Review E 2011; 84, 011305.
-  Snoeijer J.H., van Hecke M., Somfai E., and van Saarloos W. Packing geometry and statistics of force networks in granular media, Physical Review E 2004; 70:15, 011301.
-  Song C., Wang P., and Makse H.A. A phase diagram for jammed matter, Nature 2008; 453, 629–632.
-  Tatara, Y. 1989 Extensive theory of force-approach relations of elastic spheres in compression and in impact, Journal of Engineering Materials and Technology 111, 163–168.
-  Thornton C. Force Trnsmission in Granular Media, KONA Powder and Particle 1997; 15, 81–90.
-  Tighe B.P., van Eerd A.R.T., and Vlugt T.J.H. Entropy maximization in the force network ensemble for granular solids, Physical Review Letters 2008; 100:23, 238001.
-  Trappe V., Prasad V., Cipelletti L., Segre P.N., and Weitz D.A. Jamming phase diagram for attractive particles, Nature 2001; 411, 772–775.
-  Van Eerd A.R.T., Ellenbroek W.G., van Hecke M., Snoeijer J.H., and Vlugt T.J.H Tail of the contact force distribution in static granular materials, Physical Review E 2007; 75:6, 060302(R).
-  Vågberg D., Olsson P., and Teitel S. Glassiness, rigidity, and jamming of frictionless soft core disks Physical Review E 2011; 83, 031307.
-  Vicente L.N. A comparison between line searches and trust regions for nonlinear optimization, Investiga o Operacional 1996; 16, 173–179.
-  Wu, C.-Y. DEM simulations of die filling during pharmaceutical tabletting Particuology 2008; 6:6, 412–418
-  Zhou J., Long S., Wang Q., and Dinsmore A.D. Measurement of forces inside a three-dimensional pile of frictionless droplets, Science 2006; 312, 1631–1633.