Jamming transition in a two-dimensional open granular pile with rolling resistance

We present a molecular dynamics study of the jamming/unjamming transition in two-dimensional granular piles with open boundaries. The grains are modeled by viscoelastic forces, Coulomb friction and resistance to rolling. Two models for the rolling resistance interaction were assessed: one considers a constant rolling friction coefficient, and the other one a strain dependent coefficient. The piles are grown on a finite size substrate and subsequently discharged through an orifice opened at the center of the substrate. Varying the orifice width and taking the final height of the pile after the discharge as the order parameter, one can devise a transition from a jammed regime (when the grain flux is always clogged by an arch) to a catastrophic regime, in which the pile is completely destroyed by an avalanche as large as the system size. A finite size analysis shows that there is a finite orifice width associated with the threshold for the unjamming transition, no matter the model used for the microscopic interactions. As expected, the value of this threshold width increases when rolling resistance is considered, and it depends on the model used for the rolling friction.


I. Introduction
Granular materials are ubiquitous either in nature -desert dunes, beach sand, soil, etc.-or in indus-trial processes as mineral extraction and processing, or food, construction and pharmaceutical industries [1][2][3]. In fact, any particulate matter made of macroscopic solid elements can be classified as granular material. The vast phenomenology exhibited by these systems combined with an incomplete understanding about the microscopic physical mechanisms responsible for the macroscopic behavior of these materials have motivated the increasing interest of the physics community in the past years [4,5].
Although materials of this class are not sensitive to thermal fluctuations, they can be found at gas, liquid or solid phases [6]. The transition between solid and liquid phases in granular matter, which is commonly referred to as jamming/unjamming transition, has been extensively studied from both theoretical and experimental perspectives [4,[6][7][8][9][10][11]. Currently, a great effort is being made to under-stand the nature of this transition, which is still a subject of debate [12]. The jamming/unjamming transition is not a specific property of granular matter, being observed in many kinds of materials, such as foams [13], emulsions [14], colloids [15], gels [16], and also in usual molecular liquids [17] -glass transition. Liu and Nagel [9] proposed a general phase diagram as an attempt to unify the several approaches to study jamming/unjamming in disordered materials. This work has motivated several theoretical, experimental and numerical investigations, but a comprehensive understanding of this transition is still lacking.
O'Hern et al. [18,19] have performed numerical simulations of granular materials approaching to jamming in two and three dimensions. They have explored the packing fraction axis of the general phase diagram proposed by Liu and Nagel and have demonstrated, by means of finite-size analysis, the existence of a unique critical point in which the system jams in the thermodynamic limit. The authors have also shown some evidence that this point is an ordinary critical point, indicating that the jamming transition would be a second order phase transition. These results were corroborated later by experiments (c.f. Majmudar et al. [10]) and by simulations (c.f. Manna and Khakhar [20]). In Ref. [20], the authors have revealed evidence of self-organized criticality (SOC) by measuring the internal avalanches resulting from the opening of an orifice at the bottom of granular piles.
Experimental investigation of the jamming transition in granular materials under gravitational field has been conducted in a variety of ways, addressing the role of many parameters, like the grain shape, the friction coefficient, and the system geometry. However, a common feature of all these approaches is the analysis of the granular flow through bottlenecks.
The jamming of three-dimensional piles seems to be settled after the work of Zuriguel et al. [21]. They have demonstrated experimentally, for piles composed of different kinds of grains, the divergence on the mean internal avalanche size. It means that, as the outlet size approaches a critical value, the internal avalanche increases without limit and a permanent flow is established. This critical outlet is insensitive to the density, stiffness and roughness of the grains, but shows a significant dependence on the grain shape. For spherical grains, a critical out-let width w c ∼ 4.94d was obtained, for cylindrical grains w c ∼ 5.03d and for rice grains w c ∼ 6.15d.
Nevertheless, the jamming transition in twodimensional piles is still a question under debate. In order to address it, To et al. [22] have carried out experiments using two-dimensional hoppers in order to find a critical outlet size for jamming events. The jamming probability J has presented a rapid decay from J = 1 to J = 0 close to the aperture width w ∼ 3.8d, signaling a possible phase transition. The authors discussed, based on a restricted random walk model, the connection between jamming and the arch formation mechanism. Nevertheless, the point needs further investigation, especially at the limit of high hopper angle. There exist several works [23][24][25] focused on the mechanisms of arch formation, but none had explored its relation to jamming probability.
Janda et al. [26] have made some progress by simulating discharges of two-dimensional silos. They have improved the definition of jamming so that the internal avalanche size is considered. This modification addresses the extremely long relaxation times associated with jamming. Within the framework of a probabilistic theory concerning the arch formation, the authors have tested two hypothesis for the internal avalanche behavior: a functional form that predicts a divergence in the mean size, and a functional form where the mean size exists for all values of the orifice width w. Since the latter one is more compatible with the arch formation model, the authors claimed that "no critical opening size exists beyond which there is not jamming" [26].
The results mentioned so far are related to fully confined systems. Recently, a simulation study on discharges of granular piles with open boundaries [27] provided new insights on the problem. The piles are composed by homogeneous disks interacting via elastic and frictional forces. Using finite size analysis, the authors have shown that a catastrophic regime, in which the pile is completely destroyed by the opening of the orifice, is well defined. At the limit of infinitely large systems, the catastrophic regime coincides with the unjammed phase, since it implies a divergent internal avalanche. Hence, the results indicate the existence of the jamming transition. It is important to note, however, that the pile geometry could probably play a role in the causes for this distinct behavior, due to the absence of the Janssen effect, but further investigations are necessary to confirm this point.
In the present work, the investigation of jamming in 2D open systems is extended in order to consider a rolling resistance term in the grain interaction model, following the prescription adopted by Chevalier et al. [28]. The main objective of this study is the verification of rolling resistance influence on jamming. Many factors contribute to the appearance of rolling friction, including microsliding, plastic deformation, surface adhesion, grain shape, etc., but mainly, it is due to the contact deformation [29]. Here, it will be taken into account only the effect due to the contact deformation by implementing the micromechanical model proposed by Jiang et al. [30]. The rolling friction produces a resistance to roll which, among other effects, is responsible for granting more stability to granular piles [31], and for the occurrence of different types of failure modes in granular matter [32,33]. Rolling friction was also used to model a system of polygonal grains by making a correspondence between the rolling stiffness and the number of sides of the polygon [34]. It was demonstrated that it is an essential ingredient to reproduce experimental compression tests in mixtures of two-dimensional circular and rectangular grains [35]. These facts suggest that jamming could be affected by rolling friction. Nevertheless, most studies on granular materials based on computer simulations do not deal with it.
The paper is organized as follows: after a review in the Introduction, the next section is concerned with the methodology. Then, we present the results and a brief discussion. Finally, the last section gathers the conclusions and some perspectives.

II. Methods
The jamming transition is assessed by means of discharges of granular piles, simulated using the molecular dynamics method [36] with the Velocity-Verlet algorithm [37]. In a few words, the molecular dynamics consists in integrating numerically the equations of motion that governs the system dynamics. The system is constituted by N free grains governed by Newton's second law and by a finite and horizontally aligned substrate made of fixed grains. The free grains, which will form the pile, are homogeneous bi-dimensional disks that are free to translate and rotate around their center, and whose radii are uniformly distributed around an average value d, a small polydispersity of 5% was imposed in order to avoid crystallization effects. All spatial quantities will be expressed in terms of d. Since the grains have all the same mass density, their masses are proportional to their respective areas. Normalized by the mass of the heaviest grain, the masses are given by where d max stands for the diameter of the largest grain. The finite substrate of length L is composed by fixed grains of the same kind but with a smaller and fixed diameter d s = 0.1d. These grains are aligned horizontally and are equally spaced in order to form a grid without gaps.
The grains are subject to a uniform gravitational field orthogonal to the substrate line, and to short range binary interactions. Besides the viscoelastic and coulomb friction interactions used in past works [25,27,[38][39][40][41][42], two grains in contact are also subject to a rolling resistance moment due to the finite contact length l c . The rolling resistance is introduced through a micromechanical model of the contact line between two grains [30]. This model treats the contact as an object formed by a set of springs and dashpots connecting the borders of the two grains. As one grain rolls over the other, the springs in one side of the contact line contract while the springs in the opposite side stretch. This configuration generates an unbalanced force distribution and a consequent moment with respect to the grain center (see Fig. 1). This moment grows linearly as the grain rolls, until the rolling displacement δ r reaches some threshold value at which the springs located near one end of the contact line break up and new ones emerge at the other end. At that time, the moment saturates at some value that depends on the properties of the grain. The rolling displacement is defined by δ r = (ω i − ω j )∆t, where ω i refers to the angular velocity of grain i, ∆t is molecular dynamics time step, and the sum runs over time during the whole existence of contact. Based on these assumptions, the authors have derived an analytic expression for the rolling resistance moment as a function of rolling displacement. They have also proposed a simplified version -the one used in this study -in order to improve numerical computations: where k r is the rolling stiffness, µ r is the coefficient of rolling resistance, and f el is the compressive elastic force, normal to contact line. As can be noted, the expression for the rolling resistance momentum possesses a striking resemblance with the Coulomb static friction force, and as a matter of fact, the rolling resistance interaction is implemented in the molecular dynamics algorithm in the same way as the static friction. The model predicts that the rolling stiffness and the coefficient of rolling resistance are related to the contact length l c by the equations k r = k n l 2 c and µ r = µl c , where k n and µ are respectively the normal elastic constant and friction coefficient. The values used for these parameters were the same as in Ref. [27], µ = 0.5 and k n = 1000 in normalized unities (see [40] for further details). Two cases were investigated: systems in which the rolling resistance parameters k r and µ r vary according to the above-mentioned equations, and systems with fixed rolling resistance parameters, assuming that all contacts have the same deformation value.
The simulation procedure consists of two steps: (1) the formation of the granular pile with open boundaries, by deposing grains from rest, under gravity, over the substrate until an stationary state is reached; (2) the discharge itself, which consists in opening an orifice of a given width at the center of the substrate. In the first step, the initial positions of the free grains are randomly sorted along a horizontal line located at height L from the substratethe releasing height is equal to the substrate length. To avoid initial overlapping of grains, a 50% filling ratio was imposed to each line of grains released, and the time interval between successive rows is the inverse of the frequency f . Each row was released after the predecessor had fallen a distance equivalent to the maximum grain diameter. This deposition protocol mimics a dense rain of grains. During deposition, the grains may leave the system through the lateral boundaries, so that the total number of grains in the pile fluctuates as the process evolves. The release of grains ceases when the number of grains in the pile reaches a stationary value. The deposition phase ends only when a Figure 1: In panel (a), a perfectly rigid grain rolling over another one is exhibited. The contact force is concentrated on one point and is aligned through the grain center, thus, not generating momentum with respect to the center. Panel (b) shows the same situation but with deformable grains. Now, the contact force is non-uniformly distributed over a segment of length l c (contact length). The resultant force F R,c is dislocated from the center line and an opposing moment with respect to the center appears. mechanical equilibrium state is attained, and the configuration is recorded for later analysis. The equilibrium state must satisfy the following criteria: mechanical stability, absence of slipping contacts, vertical and horizontal force balance, and vanishingly small kinetic energy [43]. In the second step, the configuration recorded is loaded and an outlet of width w is opened at the center of the substrate, allowing the grains to flow through it. As the grains pass through the outlet, they are removed from the system, in order to improve computational efficiency. The simulation runs until a new equilibrium state is reached, which happens either due to the formation of an arch above the outlet or after the pile has been completely discharged. The remaining pile configuration is then recorded.
The average height h of the resulting pile after discharge, measured from the center of the substrate, was taken as an order parameter to distinguish between the two regimes. If h does not change significantly with respect to the original height, it means that an arch does readily clog the flux after the orifice opening, but if h = 0, it means that the pile has collapsed, and a jammed state was not attained for the corresponding orifice diameter and substrate length. As the orifice width w approaches a certain threshold value w t , the system suffers a transition from jammed to unjammed state. This threshold is defined as the orifice width for which the height fluctuations is maximum. The connection to the jamming transition occurs when the substrate length diverges since the collapse of an infinite substrate pile implies a continuous flowing state. Then, the jamming transition can be characterized by a critical aperture width, as defined by the following expression (2)

III. Results and Discussion
Two different models of rolling resistance interaction were tested in the simulations: the original model proposed by Jiang et al. [30] mentioned earlier, with a rolling constant and a coefficient of rolling resistance that depends on the contact length (k r = k n l 2 c and µ r = µl c ) and a simpler derived version, in which these parameters are constant over time assuming that l c = 0.05 d for all contacts. This fixed l c model assumes a mean contact length equivalent to 5% of the average grain diameter, a scenario which could be associated to a system composed by polygonal grains, for example, a extreme rolling resistance regime.
For either models, the numerical simulations described in the last section were carried out for Figure 3: Images of the pile equilibrium states for the three grain interaction models before the discharge step. The piles were grown over a L = 100d substrate and are representative samples from each type of pile. As shown in the figure, the top image represents a pile composed by grains with a fixed k r rolling resistance interaction, the image in the middle is a pile of grains with a k r ∼ l 2 c rolling resistance interaction, and the bottom image is a pile of grains without rolling resistance.
various system sizes and orifice widths. Figure  3 shows equilibrium configurations of typical piles with L = 100d for the two tested rolling resistance models and for the absence of rolling resistance case. It can be seen that the inclusion of the rolling resistance term modifies significantly the macroscopic features of the pile. It provides more stability to the structure, which is reflected by a steeper free surface, a fact also observed elsewhere [31]. Indeed, it should be noted that the rolling resistance makes possible some otherwise very unstable twograin configuration, as exemplified in Fig. 2.
The behavior of the order parameter h as function of w is presented in Fig. 4 for the three cases. Note that the transition region translates to the right as the rolling resistance becomes more important, which is an expected result since the increasing of stability allows the formation of larger arches. Figure 5 exhibits the dependence of w t (fluctuation maximum) on the system size for the three models considered and, again, the curves are dislocated along the w t axis. Nevertheless, they all share the same functional aspect, which is an evidence that the rolling resistance does not change the nature of the jamming transition, only the critical value of the threshold width. Despite the tendency to a heavy tail distribution observed for large L and w, in all scenarios, the data suggest the existence of the jamming transition. It means that the features described in Ref. [27] to characterize the transition were also observed in both scenarios tested for the rolling resistance. The fitting values were w c = (5.3 ± 0.1)d for contact length dependent k r model and w c = (7.6 ± 0.2)d for the fixed k r one, while for the model without rolling resistance, w c = (5.0 ± 0.1)d [27]. These results indicate that the rolling resistance only slightly changes the critical width when included in a system of disks, expressed by the contact length model. But, in other Fixed k r Figure 5: Graphs of the threshold orifice width w t as a function of the reciprocal system size 1/L for models with and without rolling resistance. As the caption indicates, the symbols represent the data obtained from numerical simulations and the lines are fitting curves, which provide the respective values of w c .
approaches, as in the fixed length model, it can alter significantly the critical aperture width. The probability density functions (PDF) of h for the absent rolling friction and for the fixed rolling friction models are exhibited in Fig. 6. For each model, the PDFs are obtained for three characteristic values of w: w = 1.0d , w ∼ w t , and w > w t . It can be noted that the PDFs have an approximately Gaussian peak around the initial height for all values of w, meaning that there is always a certain amount of samples that are simply not disturbed by the outlet. Apart from these samples, the great majority is completely discharged for w > w t . This result is in consonance with that obtained earlier in a different context [25]. In that study, it was shown that for large w, all blocked events occurred after only few grains had passed through the outlet. These facts indicate that the initial conditions may play an important role in jamming experiments. However, this issue needs further investigation.  the results found in granular piles without this interaction. This result strengthens the expectation that the transition exists in real granular piles and is probably affected by the system geometry. We observe that when there was rolling resistance, the piles built were more stable, denoted for the large mean height of the samples, and also more robust against perturbations, since the critical aperture width increased. In future works, we plan to present a detailed study of the arching statistics for the two approaches considered here.