Kinetically constrained model for gravity-driven granular flow and clogging

TL;DR

This paper introduces a kinetically constrained lattice-gas model with gravity to study dense granular flow and clogging, revealing how current behavior varies with density and field intensity, including jamming phenomena.

Contribution

It extends the Kob-Andersen model by incorporating gravity effects and analyzes flow and clogging behavior at different densities, providing theoretical and simulation insights.

Findings

Current matches mean field theory at low densities.

Non-monotonic current dependence on field at intermediate densities.

Flow through a narrow hole agrees with theoretical predictions.

Abstract

We add extreme driving to the Kob-Andersen kinetically-constrained lattice-gas model in order to mimic the effect of gravity on dense granular systems. For low particle densities, the current that develops in the system agrees at arbitrary field intensity with a mean field theory. At intermediate densities, spatial correlations give rise to non-monotonic dependence of the current on field intensity. At higher densities, the current ultimately vanishes at a finite, field-dependent jamming density. We supplement the study of this bulk behavior with an investigation of the current through a narrow hole. There, lateral flow decreases the local density in front of the hole. Remarkably, the current through the hole quantitatively agrees with a theoretical prediction based on the bulk current at the measured local density.

Equations

$${\eta }_{i,j}=\{\begin{array}{ll}1& \mathrm{the}\mathrm{site}\mathrm{is}\mathrm{occupied}\\ 0& \mathrm{the}\mathrm{site}\mathrm{is}\mathrm{vacant}\end{array}$$

$$J=⟨{\eta }_{0,0}(1-{\eta }_{0,1})(1-{\eta }_{-1,0}{\eta }_{0,-1}{\eta }_{1,0})(1-{\eta }_{-1,1}{\eta }_{0,2}{\eta }_{1,1})⟩.$$

$$J_{MF}(\rho )=\rho (1-\rho )(1-{\rho }^3{)}^2.$$

$$C_2(\Delta x,\Delta y)=\frac{⟨{\eta }_{i,j}{\eta }_{i+\Delta x,j+\Delta y}⟩-{\rho }^2}{\rho -{\rho }^2}.$$

$$J_{\text{tot}}=\frac{⟨N⟩}{t{L}_y},$$

$$J_{\text{tot}}=J_{\text{bulk}}(\rho ,G)w$$

$$J_{\text{tot}}=J_{\text{bulk}}({\rho }_{\text{local}},G)w.$$

Full text

I Introduction

Jamming is the transition from a flowing state to a rigid or arrested state as density increases in an assembly of particles Liu-Nature ; Hecke-JPhysCond . Colloidal suspensions, glass-forming liquids, and granular materials show heterogeneous dynamics as they approach jamming Berthier-PRE2003 ; Leonard-JCP2010 ; Pastore ; Garrahan-PRL2002 ; Berthier-JCP2003 ; Marinari-EPL ; Whitelam-PRE ; Teomy-PRE2015 . Clogging is when the flow of granular materials through a bottleneck clog-bottleneck ; Zuriguel_2005 ; Zuriguel_2011 ; Zuriguel_2012 or pipe clog-pipe cease to exist. Jamming is a bulk phenomena that occurs throughout the system, while clogging is identified as inability to flow through a specific narrow region in the system where particles get stuck. This eventually causes the flow to stop also in other regions. There have been studies of the crossover between jamming and clogging in heterogeneous environment Peter-Srep2018 . Many systems in nature and in industrial applications exhibit jamming or clogging in very complicated geometries Exp-1 ; Exp-3 , and it is important to understand how confinement induces the jamming or clogging of granular matter. Experiments on granular systems Abate-PRE show steady growth in dynamical heterogeneity as the relaxation time increases with increasing density. We employ kinetically-constrained models Sellitto_ASEP_KCM_PRL ; Shokef-EPL ; Turci_Pitard_Driving_KCM ; Teomy-PRE2014 ; Teomy-PRE2017 ; Mallick2018 to construct a lattice model for describing the physical mechanisms involved in jamming. These models do not rely on any type of interaction between the particles other than excluded volume interaction and are supplemented by specific kinetic rules that control the movement of particles depending on the local density around them. These dynamical restrictions yields glassy behavior in the system.

In this paper, we examine the effect of an external field on jamming in the two-dimensional Kob-Andersen Kob-PRE lattice gas with kinetically-constrained dynamical rules. We implement the external field in a way that, to our understanding, better represents gravity-driven flow of granular materials, compared to previous studies Sellitto_ASEP_KCM_PRL ; Turci_Pitard_Driving_KCM . We measure the bulk current as a function of bulk density and observe that the current vanishes at a jamming density which depends on the applied external field. We also measure the bulk current as a function of the external field for the whole range of bulk densities. To understand the microscopic origin for the dynamically jammed states, we find that these states exhibit some rare mobility regions along with other regions of higher mobility due to local relaxation. Increasing or decreasing local mobility in some region sometimes facilitates the overall cooperative dynamics of the system. Furthermore, we study the effect of confinement on clogging by investigating the behavior of the system when particles are allowed to pass through a narrow orifice. We quantitatively explain the measured current in terms of the current of an effective bulk system with the local density that we measure in front of the orifice.

The paper is structured as follows. In Sec. II, we introduce our model and describe the numerical simulation details. In Sec. III, we study the bulk current as a function of external field and show the non-monotonic behavior of current. We measure two-point density correlations, focus on the microscopic structure of the system at different parameter regime, and obtain the jamming phase-diagram vs density and field. In Sec. IV, we study our model with confinement, which mimics the flow of particles through an orifice. Section V summarizes the work.

II Model

We use the Kob-Anderson kinetically-constrained model Anderson-PRL , in which particles are randomly distributed on the square lattice with density $\rho$ and they randomly attempt to move in one of four possible directions. According to this model, the move is allowed only if the target site is vacant and if there are at least two vacant nearest neighbors before and after the move. An external field may be applied in some preferred direction in a manner, which we will refer to as the E-Model Sellitto_ASEP_KCM_PRL ; Turci_Pitard_Driving_KCM . Then, each particle has a rate $R$ to move along the field or in each of the transverse directions, whereas there is a reduced rate $Re^{-E}$ for particles to move against the external field, see Fig. 1a.

We suggest that in gravity driven granular flow, the rate of motion along the field should be substantially higher than in the transverse directions. Thus we introduce the G-model with a constant rate $R$ to move along the field together with transverse movement with rate $Re^{-G}$, and no movement against the field, as shown in Fig. 1b. Here particles have a constant rate to attempt moving in the direction of the field $G$, and particle current changes by varying the lateral motion, which in turn influences the acceptance rate of moves along the field. Note that the extreme driving limit of the E-model, $E=\infty$, is equivalent to $G=0$ in the G-model. Thus by increasing $G$, we are changing the rate of motion in the transverse direction and thus extending the E-model to even more extreme driving.

III Bulk Current

In order to understand the behavior of the G-model precisely, we study the dynamics of the system in various regimes. Here we begin with the bulk behavior and consider a system with periodic boundary condition in both directions. In the low-density regime, we expect correlations to be weak and therefore we can predict the current of the particles in the direction of the field. We define the occupancy index $\eta_{i,j}$ for each site $(i,j)$ in the two-dimensional lattice as:

[TABLE]

Now we can write the mean current in the direction of the field, based on the kinetic constraints. It can be clearly seen from Fig. 2 that when a particle at $(0,0)$ attempts to move to the vacant site $(0,1)$, then at least one of the sites denoted by light blue must be empty along with at least one of the sites denoted by light green that must also be empty to satisfy the kinetic constraint. Therefore the current reads,

[TABLE]

Without correlation we obtain the mean-field (MF) current,

[TABLE]

We run dynamical Monte Carlo simulations with a rejection-free algorithm, which makes the simulations very efficient, on a periodic lattice of dimension $L\times L$ and typically use $L=400.$ We first allow the system to relax for $t=10^{6}$ time steps and then start measuring the current after the system has reached the steady state. We compare our numerical result with the MF current, first in the two extreme limits: $G=0$ and $G=\infty$, having only longitudinal motion in the later case.

Figure 3 shows good agreement with the MF approximation, given in Eq. (2) in the low-density regime. As follows from this MF approximation, there is no dependence on the field $G$ at all. Above $\rho\approx 0.2$, correlation start playing a substantial role and the observed current is lower than the MF prediction. As $G$ increases, the rate to move in the direction transverse to the field is lower, thus creating more constrained dynamics, leading to stronger deviation from MF behavior. We also observe that, although our MF theory predicts $J>0$ for any density $\rho<1$, in our simulations the current stops at some finite density, which is dependent on the field $G$. Specifically, for $G=0$, jamming occurs at $\rho=0.83$ and for $G=\infty$, it occurs at $\rho=0.78$. This jamming at finite density has been observed in the $E=\infty$ limit of the E-model Sellitto_ASEP_KCM_PRL , which coincides with our $G=0$ limit.

Note that these values of the jamming density are significantly lower than what we would expect in the Kob-Andersen model at such system sizes without driving. Although in the thermodynamic limit the undriven model jams only at $\rho=1$ TBF , finite size effects lead to jamming at $\rho_{c}=1-\frac{\lambda}{\log L}$, where $\lambda$ may be approximated by $\lambda=0.25$ Teomy_JCP . Thus for our $L=400$ system size, this would lead to $\rho_{c}=0.96$.

III.1 Current as a function of $G$

We simulated the steady state current behavior as a function of the external field $G$. The results of the simulation for varying $G$ values can be seen in Fig. 4. For each density, we normalize the current by the MF prediction according to Eq. (2). As we have already seen, the substantial deviation from the MF result, and the significant dependence on $G$, both take place at higher densities. This can be understood by the fact that the rate to move sideways controls the probability to get stuck in the way. Whenever the rate to move sideways is large ($G$ is low), there is a higher probability that a particle will change it’s trajectory and hence will be less correlated, causing the current to increase at lower values of $G$.

Remarkably, in a finite region of densities, $0.66\leqslant\rho\leqslant 0.76$, the dependence on $G$ becomes non-monotonic and there is a region where the current slightly increases as we increase $G$ until the current reaches its asymptotic value at infinitely large values of $G$.

III.2 Jamming phase diagram

According to out MF prediction, the current asymptotically vanishes only as the density approaches unity. However, from the numerical observations, one can see that there is some field-dependent critical density $\rho_{c}(G)$, at which the particles are arrested and the system gets jammed. Specifically, Fig. 3 shows that the current goes to zero at different densities for $G=0$ and $G=\infty$. This implies that such jamming effect depends on the field intensity $G$. In order to investigate the transition from free flow to a jammed state, we obtained the phase diagram shown in Fig. 5, which describes the numerical results of the current as a function of the density and field.

As we can see from the phase diagram, the transition to the jammed state happens at the density range $0.78<\rho_{c}(G)<0.83$. This critical density was observed in the past for the $E=\infty$ limit of the E-Model Turci_Pitard_Driving_KCM , which is equivalent to $G=0$ in our G-Model. However, we have achieved the extension to higher values of $G.$ We find that as the density increases, a smaller field is required to get the jammed phase.

III.3 Correlation analysis

We have performed correlation analysis by defining the two-point density correlation function as,

[TABLE]

As can be seen from Fig. 6, the correlations are mostly longitudinal as $G$ becomes stronger. This is because increasing $G$ is equivalent to decreasing the probability to move in the transverse direction, which corresponds to the dominant motion being longitudinal.

III.4 Analysis of microscopic structure

One of the prominent differences between various regimes of $G$ values are the internal structures that are formed in the steady state. As the density increases, the homogeneously distributed particles start to create clusters of particles and clusters of holes, as can be observed from Fig. 7. At the high densities, the typical shapes of the microscopic structures of vacancies are illustrated in Fig. 8b. The formation of these hole clusters starts with the formation of shear bands Turci_Pitard_Driving_KCM which are perpendicular to the applied field. These shear bands basically act as barriers to the flow, below which there are areas of vacancies, which go around these barriers and, due to transverse motion, gradually fill in these areas of vacancies.

As can be seen from Fig. 8a, for small values of $G$, the height of the vacancies region is typically small, as there is a high transverse motion, and the particles fill the open space below the shear bands relatively quickly, as opposed to the typical structure of the hole cluster for high values of $G$, as illustrated in Fig. 8b. As we further observe, these hole clusters are not persistent and they dynamically form and get destructed in the steady state. The evolution of the configuration at the steady state and the appearance of the hole clusters can be seen in Fig. 7. This effect starts taking place at high densities. For $G=\infty$, there is no motion at all in the transverse direction and the only motion is in the direction of the field. In this case, it can be seen that the hole clusters are the longest, as depicted in Fig. 8b.

IV Flow through an orifice

In this section, we study our model with confinement that mimics the behavior of discharging granular materials through an opening. The confinement reveals new phenomena, such as density heterogeneity and clogging of the system near the opening. In addition, we would like to investigate the confined flow and to reveal its connection to the bulk current, discussed in Sec. III.

We define a confined geometry imitating the discharge of grains through an orifice, as shown in Fig. 9. Also in this geometry, we employ periodic boundary condition in both directions, whereas, at the center of the lattice we place a rigid horizontal barrier, with a narrow opening of width $w$. The obstacle is achieved by artificially placing one row of occupied sites in the initial condition and not moving them during the simulation. Due to the periodic boundary conditions, this geometry is equivalent to infinitely many openings with distances $L_{x}-w$ in the transverse direction, and infinitely many such barriers at a distance $L_{y}$ in the longitudinal direction. In order to eliminate the finite size effect, we maintained a considerable value of the ratio $L_{x}/w$.

The simulation begins with randomly and uniformly distributed particles throughout the lattice. Due to the presence of the barrier, there is sedimentation of the particles towards the barrier, whereas in front of the opening there is flow. As a result, heterogeneity of the density develops. For $G<\infty$, due to lateral motion, particles rearrange, such that the density above the barrier increases, and by mass conservation, the density above the opening decreases, as shown in Fig. 10. We also observe that in the steady state, the only contributing part to the current is the flow above the orifice, denoted by the yellow particles in Fig. 10. This is because the system is jammed in other regions of the lattice, as shown by the light blue particles.

We are interested in measuring the total current of the particles through the opening, namely the number of particles flowing through the opening per unit time as given by,

[TABLE]

where $\left\langle N\right\rangle$ is the average number of particles that moved along the direction of the field during the time interval $t$ anywhere in the lattice, and $L_{y}$ is the lattice height (see Fig. 9). Using this definition, we can now estimate the total current through the orifice using the bulk current, which we studied in Sec. III. Recalling that the bulk current is defined as the number of particles flowing in the system per unit time and per lattice size, we can formally define the bulk current as, $J_{\textrm{bulk}}=\left\langle N\right\rangle/(tL_{x}L_{y}).$ In order to be consistent with the definition of the total current in Eq. (4), our crudest approximation will be as follows,

[TABLE]

Figure 11 shows results of numerical simulation, where we calculated the total current through the orifice, in comparison to the crude approximation given in Eq. (5). As we can see from the plot, there is a considerable discrepancy between the numerical results (blue points) and this calculation (red line). To explain this, we recall that the density in the confined system is heterogeneous, and in particular, the density above the orifice is lower than the average density in the system. Hence we consider the local density measured above the orifice, as this is the only region that contributes to the current in the steady state. This corrects our theoretical prediction to,

[TABLE]

In Sec. III, we have seen that in non-confined systems, the bulk current is higher for lower densities. Remarkably, the analytical expression (6) that takes into account the reduced local density above the opening for the total current coincides with the numerical results, as seen in Fig. 11. Thus, using the behavior studied for the bulk system, we can predict the dynamics in confined systems. We can perform a similar comparison also for different values of $G.$ As the field increases, the transverse motion becomes less probable due to a lower rate to move in that direction. This leads to a longer time needed to reach the steady state in confined systems. In Fig. 12 we show the numerical results for two different values of $G$ with remarkable agreement with the same analytical prediction, described by Eq. (6). As $w$ becomes smaller, the current decreases linearly with $w$, as Eq. (6) predicts. For very small values of $w$, we obtain clogging of the particles near the opening and the overall current ceases for some realizations. The numerical results shown in Fig. 11 and 12 are obtained by averaging over the current from many realization that weren’t clogged. This clogging state becomes more and more probable as $w$ becomes smaller, and it would be interesting to study it further.

V Conclusions

We studied the behavior of granular materials subjected to a gravitational field in a confined geometry. To that end, we used the kinetically-constrained Kob-Anderson model, which implements simple dynamics of a lattice gas, with no interactions between neighboring particles, but with kinetic rules that depend on the occupancy of neighboring sites. Such a description allows us to employ efficient numerical simulations. In order to simulate the gravitational field, we modified the rates of the moves in the direction perpendicular to the field, controlled by the parameter $G$, and did not allow moves against the direction of the field. We first studied the bulk behavior of the model, eliminating confinement by creating a large system with periodic boundary conditions. We built a MF theory and performed numerical simulations in order to understand the range of validity of the MF theory and the role of correlations. For that purpose, we constructed an efficient rejection-free algorithm, in which we skip all the unsuccessful moves. We observed that at the range of densities $0.66\leqslant\rho\leqslant 0.76$ the current exhibits non-monotonic behavior as a function of the driving field $G$. From the microscopic point of view we see that at the high density regime, close to jamming, the spatial distribution of the particles is very heterogeneous, namely there are spatial structures of holes that are dynamically formed and destructed as a result of the development of shear bands.

To study discharge through a narrow opening, we defined a confined geometry in our 2D lattice by introducing a horizontal barrier with an opening of some size that we could change. We conducted numerical simulations, measuring the total current through the opening, and revealed two phenomena - clogging close to the opening and highly reduced density of the flowing particles. We measured the spatial distribution of the particles in the steady state with confinement, measuring only the current above the opening, as this is the only contributing part to the current in the whole system. We were able to analytically predict the current as a function of the opening width by using the bulk current corresponding to the reduced local density, as measured from numerical simulations with confinement. We obtained excellent agreement between the total current through the opening from numerical simulation and the analytical calculation based on the current measured in the bulk case, corresponding to the reduced local density above the opening.

One of our most remarkable observations for bulk behavior is the non-monotonicity of the current as a function of field intensity at intermediate densities. Further structural analysis could identify the reason for this behavior.

Similar studies of discharge through an opening and related phenomena such as spatial density heterogeneity and bottleneck clogging have already been seen in granular matter, whether in experiments or simulations. Our model reveals qualitatively the same phenomena, whereas thanks to its simplicity the simulations are fast and can be done on very large systems. The model can also be extended to three dimensions where, instead of looking on the square lattice, one may consider the cubic lattice. This case may perhaps be more suitable for comparison with actual granular experiments.

VI Acknowledgments

We thank David Gomez, Erdal Oguz, Daniel Vilyatser, Eyal Rubin, Nimrod Segall, Hadas Shem-Tov, and Eial Teomy for helpful discussions. This research was supported by the Israel Science Foundation Grant No. 968/16 and by the Prof. A. Pazy Research Foundation.