Structural Evolution of Granular Systems: Theory

TL;DR

This paper develops a comprehensive theoretical framework for understanding how the structure of planar granular systems evolves over time, including steady states and stability, validated by simulations.

Contribution

It introduces a general theory for cell order distribution evolution in granular systems, including dynamic equations and stability analysis, extending to structural organization.

Findings

All steady states are stable.

Steady states satisfy detailed balance for cell order ≤ 6.

The theory is validated against extensive simulations.

Abstract

A general theory is developed for the evolution of the cell order (CO) distribution in planar granular systems. Dynamic equations are constructed and solved in closed form for several examples: systems under compression; dilation of very dense systems; and the general approach to steady state. We find that all the steady states are stable and that they satisfy detailed balance-like condition when the CO$\,\leq 6$. Illustrative numerical solutions of the evolution are shown. Our theoretical results are validated against an extensive simulation of a sheared system. The formalism can be readily extended to other structural characteristics, paving the way to a general theory of structural organisation of granular systems.

Figures

Equations

$$\dot{n_k}$$

$${\eta }_{i,j}\equiv p_{i,j}{Q}_i{Q}_j-q_{i+j-2,i}{Q}_{i+j-2}.$$

$$\dot{Q_k}$$

$$-\sum _{i=k+1}^{\infty }{\eta }_{k,i-k+2}\left(1+{\delta }_{i,2k-2}\right)+Q_k\sum _{\genfrac{}{}{0ex}{}{\text{all possible}}{\text{processes i,j}}}{\eta }_{i,j},$$

$${\dot{N}}_c=\sum _{k=3}^{\infty }{\dot{n}}_k=-N_c\sum _{\genfrac{}{}{0ex}{}{\text{all possible}}{\text{processes i,j}}}{\eta }_{i,j}.$$

$$\frac{{\dot{N}}_c}{N_c}=\frac{\dot{z}}{z-2}.$$

$${\dot{Q}}_3$$

$${\dot{Q}}_4$$

$$\dot{z}$$

$$t-t_0=\frac{1}{p_{3,3}}\ln \left[{\left(Q_3-a\right)}^{\alpha }{\left(Q_3-b\right)}^{\beta }{\left(Q_3-2\right)}^{\gamma }\right],$$

$$\delta \dot{\vec{Q}}(t)=A\cdot \vec{Q}(t),$$

$$\vec{Q}(t)={\vec{Q}}^s+\sum _i({\vec{v}}_i^{†}\cdot \delta \vec{Q}(t_0))e^{{\lambda }_{i}t}\cdot {\vec{v}}_i$$

$$(Q_3^s{)}^3+{\theta }_{3,4}(Q_3^s{)}^2+{\theta }_{3,3}{\theta }_{3,4}{Q}_3^s={\theta }_{3,3}{\theta }_{3,4},$$

$$Q_4^s=\frac{1}{{\theta }_{3,3}}{Q}_3^s;Q_5^s=1-Q_3^s-Q_4^s.$$

$$\delta {\dot{Q}}_k$$

$$-\sum _{i=k+1}^{i_{\mathrm{m}}}\left[p_{k,i-k+2}(Q_k^s\delta Q_k+Q_{i-k+2}^s\delta Q_{i+k-2})-q_{i,k}\delta Q_i\right](1+{\delta }_{i,k-i+2})$$

$$-Q_k^s\sum _{\ell =3}^{\mathcal{C}}$$

$$-\sum _{i=\ell +1}^{i_{\mathrm{m}}}\left[p_{\ell ,i-\ell +2}(Q_{\ell }^s\delta Q_{\ell }+Q_{i-\ell +2}^s\delta Q_{i+\ell -2})-q_{i,\ell }\delta Q_i\right](1+{\delta }_{i,\ell -i+2})\},$$

$$\sum _{k=3}^{\mathcal{C}}\delta {\dot{Q}}_k$$

$$-\sum _{i=k+1}^{i_{\mathrm{m}}}\left[p_{k,i-k+2}(Q_k^s\delta Q_k+Q_{i-k+2}^s\delta Q_{i+k-2})-q_{i,k}\delta Q_i\right](1+{\delta }_{i,k-i+2})\}$$

$$-\underset{=1}{\underbrace{\sum _k{Q}_k^s}}\sum _{\ell =3}^{\mathcal{C}}$$

$$-\sum _{i=\ell +1}^{i_{\mathrm{m}}}\left[p_{\ell ,i-\ell +2}(Q_{\ell }^s\delta Q_{\ell }+Q_{i-\ell +2}^s\delta Q_{i+\ell -2})-q_{i,\ell }\delta Q_i\right](1+{\delta }_{i,\ell -i+2})\}.$$

$$\sum _{k=3}^{\mathcal{C}}\delta {\dot{Q}}_k$$

$$-\sum _{i=k+1}^{i_{\mathrm{m}}}\left[p_{k,i-k+2}(Q_k^s\delta Q_k+Q_{i-k+2}^s\delta Q_{i+k-2})-q_{i,k}\delta Q_i\right](1+{\delta }_{i,k-i+2})\}$$

$$-\sum _{k=3}^{\mathcal{C}}$$

$$-\sum _{i=k+1}^{i_{\mathrm{m}}}\left[p_{k,i-k+2}(Q_k^s\delta Q_k+Q_{i-k+2}^s\delta Q_{i+k-2})-q_{i,k}\delta Q_i\right](1+{\delta }_{i,k-i+2})\}$$

$$=0\text{.}$$

Full text

Granular dynamics, structural evolution, cell order distribution

Introduction: Understanding and modelling self-organisation of dense granular matter (DGM) under external forces is essential to many natural phenomena and technological applications. Examples are: consolidation and failure of soils, packing of particulates in technological processes, initiation of avalanches, flow of slurries, and dense colloidal suspensions, to mention a few. This is also one of the most important problems in granular science EdOa89a ; Vogel_Roth_2003 ; Cheng_etal_1999 ; BaBl02 ; Aste_etal_2007 ; Meyer-etal2010 because both the dense flow dynamics and the large-scale properties of consolidated DGM depend strongly on the particle-scale structure. Such properties are: permeability to fluid flow Vogel_Roth_2003 , catalysis, heat exchange, and functionality of fuel cell electrodes Ametal17 . In this paper we address this generic problem and develop a set of equations for the quasi-static structural evolution of two-dimensional (2D) systems of rigid grains.

Structural evolution of DGM is mediated by continual making and breaking of intergranular contacts. In turn, the temporal structural configuration determines the force transmission in the medium, which drives the evolution. The intergranular contact network can be regarded as a graph containing voids (or cells), each of which characterized by the number of grains enclosing it – its order. Developing a theory for the evolution of the cell order distribution (COD) is the specific aim here. This distribution has been argued Fretal08 and shown MaBl14 ; MaBl17 to converge to a universal form, in which the intergranular friction is scaled away. Specifically, a set of evolution equations is constructed and solved in closed form, under some assumptions, both for very dense closed systems and for closed systems approaching a limit steady state. More general cases are solved numerically. We then compare the theoretical results with numerical simulations of sheared systems and show that there is good agreement.

The evolution equations: When a 2D DGM undergoes slow quasi-static dynamics, it is possible to identify, at every moment, a contact network. In this graph the grain centres (however defined) are the vertices and edges connect grains in contact, with each edge corresponding to a specific contact. These make cells and we assign each a cell order defined as the number of grains surrounding it. In such dynamics, the stress response is faster than any other process and a static stress state is established practically immediately after a contact is made or broken. The systems can be regarded then as moving from one stress state to another.

We consider systems free of gravity and body forces. In the concluding discussion, we argue that including these simplifies the following formalism. The structure evolves through making and breaking of these contacts, which we call contact events (CEs). The former splits a cell into two smaller ones and the latter merges two cells into a larger one. If a CE involves two rattler-free neighbour cells of orders $i$ and $j$, the two processes satisfy $i+j\xrightleftharpoons{}i+j-2$, exemplified in Fig. 1. Analogously, if the mother cell contains a rattler that participates in the event, the process sum rule is $i+j\xrightleftharpoons{}i+j-4$. To model the dynamics of the COD, we define $n_{k}$ as the total number of $k$-cells in the system and their fraction of the total cell population, $N_{c}$, as $Q_{k}=n_{k}/N_{c}$. We also define the following rates, which we assume are system size-independent: $p_{i,j}$ = the merging rate of an $i$- and a $j$-cells into an $i+j-2$-cell; $q_{k,i}$ = the splitting rate of a $k=i+j-2$-cell into an $i$- and a $k-i+2$-cells; $r_{i,j}$ = the merging rate of an $i$- and a $j$-cells, containing a rattler, into a $k=i+j-4$-cell; and $s_{k,i}$ = the splitting rate of a $k$-cell, containing a rattler, into an $i$- and a $(k-i+4)$-cell. The COD evolves via four basic CEs. Two of creation: a $k$-cell is either a merger of $j$ and $k-j+2$, or an offspring of a split large-order cell, and two of annihilation: either by splitting into two offsprings or merging to make a larger cell. Each of these processes has an equivalent when the combined cell contains a rattler, in which case the rates $p_{j,k-j+2}$ and $q_{k,j}$ are replaced, respectively, by the rates $r_{j,k-j+4}$ and $s_{k,j}$.

If very large cells are rare then the occurrence of cells containing more than one rattler is very low and we ignore such events in the following. Then the evolution equations are:

[TABLE]

The 1/2 factor corrects for double counting and the terms containing $\delta$-functions account for loss (gain) of two same order cells upon creation (annihilation) of a larger cell. The ’tilded’ parameters in (1) are related directly to the size-independent rates: $\tilde{p}_{i,j}\equiv p_{i,j}/N_{c}$, $\tilde{q}_{k,i}\equiv q_{k,i}/N_{c}$, $\tilde{r}_{i,j}\equiv r_{i,j}/N_{c}$, $\tilde{s}_{k,i}\equiv s_{k,i}/N_{c}$. To simplify the following analysis, we ignore rattlers and set $r_{i,j}=s_{i,j}=0$. Including the more realistic rattler-related terms is straightforward, but it would result in cumbersome expressions without adding any more insight. This amounts to assuming that the average number of rattlers per cell type is constant and therefore so is the number of non-rattlers. This assumption is indeed borne out by our numerical simulation results.

From eqs. (1), we note that $\sum_{k}(k-2)n_{k}=E$ is a conserved quantity. This quantity has a physical interpretation: $\sum_{k}n_{k}\equiv N_{c}$ and $\sum_{k}kn_{k}\equiv Nz$ is twice the number of contacts, with $z$ as the mean coordination number and $N$ as the number of grains. We use Euler’s topological expression in the plane to relate the numbers of vertices (=$N$ grains), edges (=$Nz/2$), and cells, $N_{c}$, $N-Nz/2+N_{c}=\mathcal{O}(\sqrt{N})$, in which $\mathcal{O}(\sqrt{N})$ are boundary terms. It follows that $E=2N+\mathcal{O}(\sqrt{N})$. Embedding the system on the surface of a sphere, $E=2(N-1)$ exactly.

To make (1) size-independent, we convert them for the fractions $Q_{k}$. To this end, we define the deviation of the rattler-free steady state process $i+j\xrightleftharpoons{}i+j-2$ from ‘detailed balance’-like steady state,

[TABLE]

When $\eta_{i,j}=0$, the ‘reaction’ $i+j\to i+j-2$ and ‘back-reaction’ $i+j-2\to i+j$ occur at the same rate, satisfying the definitions of Lewis and Ter Haar for detailed balance in thermodynamic systems DBL . The sign of $\eta_{i,j}$ determines which direction of the process $i+j\leftrightharpoons i+j-2$ is more frequent. Generically, dilation or compression correspond to positive and negative $\eta_{i,j}$, respectively. Similarly, $\zeta_{i,j}\equiv r_{i,j}Q_{i}Q_{j}-s_{i+j-4,i}Q_{i+j-4}$ is the deviation of the rattler-involving steady state process $i+j\xrightleftharpoons{}i+j-4$ from a detailed balance-like state.

Noting that $\dot{Q_{k}}=\left(\dot{n_{k}}-Q_{k}\dot{N_{c}}\right)$, the equations for the cell fractions become

[TABLE]

in which we used

[TABLE]

Eqs. (3) are now conveniently system size-independent. In practice, cell orders in realistic quasi-static systems cannot exceed an upper bound $\mathcal{C}$ and remain mechanically stable.

Focusing on closed systems with $N(\gg 1)$ constant, we use Euler’s relation again to derive a relation between the rates of change of $N_{c}$ and $z$:

[TABLE]

Exact solutions for $3$-$4$ systems: The simplest non-trivial systems to study, realisable under high-compression processes, consist only of $3$- and $4$-cells. Eqs. (3) then reduce to

[TABLE]

Since $Q_{3}+Q_{4}=1$ the first two equations are dependent. Integrating (6) and using (2) yields (see supplementary material)

[TABLE]

with $a$ and $b$ the roots of $p_{3,3}Q_{3}^{2}+q_{4,3}Q_{3}-q_{4,3}$, and $\alpha=[(a-2)(a-b)]^{-1}$, $\beta=[(b-2)(b-a)]^{-1}$, and $\gamma=[(a-2)(b-2)]^{-1}$. From this solution and (6) we can obtain $Q_{4}$ and $z$. Examples of these solutions are shown in Fig. 2. We observe that the steady state is the same regardless of the different initial states and is determined only by the rates. As we show explicitly below, for $\mathcal{C}\leq 6$, this is not only a general feature but the steady state is also uniquely determined by the rates. This is illustrated for several $3$-$4$-$5$ systems in Fig. 3 (b) and discussed in more detail below.

The steady state: The rightmost sum in (3) describes the rate of change of the total number of cells and it vanishes at steady state. A direct calculation of eqs. (3) for dense systems provides two significant observations: for systems consisting of highest cell order, $\mathcal{C}\leq 6$., (i) the steady state is unique; (ii) this steady state satisfies $\eta_{i,j}=0$ for all $i,j$, which implies a detailed balance-like state – a very surprising result in systems that are manifestly far from conventional equilibrium.

For systems with $\mathcal{C}>6$, there are infinitely many steady state solutions, in addition to the detailed balance-like one (see details in the supplementary). For example, steady states of systems with highest order $\mathcal{C}=7$ need only satisfy $\eta_{4,5}=\eta_{3,6}=0$ and $\eta_{4,4}=-\eta_{3,4}=-\eta_{3,5}=\eta_{3,3}$. Therefore, if $\eta_{3,3}\neq 0$, such steady states exist and they are not detailed balance-like.

The approach to steady state: It is useful to understand the approach to the steady state, as well as its stability. Near this state, $Q_{k}(t)=Q_{k}^{s}+\delta Q_{k}(t)$ is only slightly different from its steady state value, $Q_{k}^{s}$, with $\lvert\delta Q_{k}\rvert\ll Q_{k}^{s}$. Defining $\vec{Q}(t)=\vec{Q}^{s}+\delta\vec{Q}(t)$, with $\vec{Q}\equiv(Q_{3},\dots,Q_{\mathcal{C}})$, and expanding the r.h.s of eq. (3) to linear order, we obtain for $\delta\vec{Q}(t)$

[TABLE]

in which the components of the constant matrix $A$ are cumbersome combinations of $p_{i,j}$, $q_{i,j}$, and the steady state fractions $Q_{k}^{s}$.

Denoting the $i$th eigenvector and eigenvalue of $A$, respectively, by $\vec{v}_{i}$ and $\lambda_{i}$, we have near the steady state

[TABLE]

with some initial time $t_{0}$. As we show in the supplementary material, all $\lambda_{i}$ must be real, there are no strictly positive eigenvalues and there has to be at least one negative eigenvalue.

To determine the unique steady-state solution of the ${3\text{-}4}$ system, we use the normalisation and detailed balance-like, $\eta_{3,3}^{s}=0$, condition to find ${Q_{3}^{s}=[-\theta_{3,3}+\sqrt{\theta_{3,3}^{2}+4\theta_{3,3}}]/2}$, with $\theta_{i,j}\equiv q_{i+j-2,i}/p_{i,j}$. The normalisation confines the dynamics to a line in the $Q_{3}$-$Q_{4}$ plane, with the steady state as unique stable fixed point on it that is independent of the initial state. This was tested numerically for several systems and is illustrated in Fig. 3(a), in which the arrow lengths represent the rate of approach to the steady state for $p_{3,3}=q_{4,3}=1$.

The dynamics of $3$-$4$-$5$ systems can be analysed similarly. Using $\eta_{3,3}^{s}=\eta_{3,4}^{s}=0$ and the normalisation condition, the steady state solution satisfies

[TABLE]

and

[TABLE]

Using the positivity of the $Q_{k}$, the uniqueness of the steady state can be established by studying the position of extrema of (10) WaThesis18_253 and it is determined by the rate fractions, $\theta_{i,j}$. The normalisation confines the dynamics to the plane $Q_{3}+Q_{4}+Q_{5}=1$ and the approach to steady state in this plane is illustrated in Fig. 3(b) for two different sets of rates, which give two distinct steady states.

For $\mathcal{C}\geq 6$, the dynamics take place on the hypersurface $\sum_{k}Q_{k}=1$, on which there is at least one detailed balance-like steady state. As discussed above, for $\mathcal{C}>6$, there may be additional steady states. Moreover, since all the steady states are stable, then, if they exist, they must be distinctly separate on the hypersurface.

Simulations: To test the theory, we carried out numerical simulations of quasi-static two-dimensional simple shear between parallel plates. We use a dissipative spheres model Bretal96 ; Sietal01 and the implementation of the model was provided by the open source project LIGGGHTS Kletal12 . Within this implementation, we used a Hertz force model for the grain-grain and grain-walls interactions, with: Young’s Modulus $5\times 10^{6}$Pa, restitution coefficient $0.2$, Poisson’s ratio $0.5$, and friction coefficient $0.5$. The system consisted of $N=21,690$ spheres restricted to move in the plane $z=0$ of four different diameters, $7$mm ($5402$ grains), $9$mm ($5400$), $11$mm ($6120$) and $14$mm ($4770$), respectively.

The system has a length of $1.65$m and periodic boundary conditions in the $x$-direction. The initial state was generated by compressing the grains from an initially loose random distribution by a flat surface at constant pressure, $P_{\mathrm{initial}}=54.45$Kg/m, in the $y$-direction until their total kinetic energy fell below $10^{-12}$J per particle. This state was regarded as mechanically stable. Then, maintaining the confining pressure, shear started in the $x$ direction, with shear velocity $v_{\gamma}=0.06$m/s. No gravitational field was applied.

The time-step was $10^{-6}$s and the grain positions and velocities were saved every $200$ time-steps for collecting detailed data on the contact network evolution.

At each stop, cells and the grains surrounding them were identified, from which we tracked the evolution of the COD and obtained the rates of contact events, $p_{ij}$ and $q_{ij}$. Cells of orders higher than $13$ occurred very rarely and, therefore, were excluded from the analysis.

Analysis of simulation data: The rates are determined from the numbers, $N^{(p)}_{i,j}\equiv p_{i,j}\,Q_{i}Q_{j}N_{c}\Delta t$ and $N^{(q)}_{i,j}\equiv q_{i+j-2,i}\,Q_{i+j-2}N_{c}\Delta t$, of the events ${i+j\xleftrightharpoons{}i+j-2}$ during a time interval $\Delta t$. In general, the rates should depend on the force distribution and the shear rate, both of which change slightly during the simulation as the system dilates before reaching a steady state. Since the above analysis is for constant rates, this could complicate a direct comparison between the analytical solutions and the simulation data. Fortunately, the rates change little and their time dependence can be neglected. Fewer than $15$ cells of a specific CO or fewer than $40$ events per $0.1$s of a specific process $\eta_{i},j$ lead to large statistical errors and were ignored. Interestingly, in spite of the long simulation time, the rates fluctuated significantly, $30\%-65\%$, perhaps because of the sensitively to small changes in $Q_{k}$. Smaller fluctuations may be achieved in larger systems with more contact events. Fig. 4 shows the rate diagram, calculated from data collected during intervals of $0.1$s and averaged over the entire evolution process. Such rate diagrams characterise uniquely a dynamic process.

Using the calculated rates, shown in Fig. 4, we solved the evolution equations (3) numerically, with the initial COD determined by averaging over the first $0.1$s in the simulation. While there was a reasonably good agreement between the solution and the simulation data when using the computed mean rates, we found that an even better agreement was achieved by correcting only $q_{6,4}$ by $10\%$. Note that this correction is much smaller than the aforementioned statistical fluctuations of the rates. The good agreement between our solution and the simulation data after this minor correction is shown in Fig. 5.

Conclusion: To conclude, we have constructed master equations to describe the evolution of a key structural characteristic of granular media - the cell order distribution (COD) - from any initial state, given the cell merging and splitting rates. The equations yield a surprising result: the steady states of the non-equilibrium dynamics of granular systems, with cell order no higher than $\mathcal{C}=6$ satisfy a detailed balance-like condition. This detailed balance-like steady state is unique and stable. Systems including cell orders of $7$ and higher can converge to this steady state, as well as to other solutions that do not support detailed balance, all of which are stable fixed points of the equations.

We validated the theory by running a long simulation of a sheared system, determining the rates, and using those to solve the master equations. Indeed, the calculated solution agrees nicely with the simulated COD evolution, in spite of large fluctuations in the rates computed from the simulations. These fluctuations may be caused by ‘clappers’: Particle pairs that make and break their common contact repeatedly. We plan an experiment to study the relative contributions of clappers and non-clappers to the rate statistics.

In our analysis, we assumed constant rates throughout the dynamic process. This is unlikely to be the general case and extending the theory to time-dependent rates is the next step. Such time dependence is expected because the contact event rates is likely to be sensitive to the intergranular force distribution, which is position- and time-dependent. Indeed, it has been argued Bletal15 ; MaBl17 that the structural evolution of granular system is a self-organisation process coupled to the evolution of the intergranular forces. Therefore, these equations and their extension form an important step towards a complete understanding of the structure-forces co-evolution. Another extension would be to open systems, when the number of particles is not conserved. This extension could be related to a grand-canonical statistical mechanical description of granular ensembles Bl16 . Although the evolution eqs. (1) include effects of rattlers, we disregarded these in the analysis for clarity. Indeed, very dense systems with low values of $\mathcal{C}$ should be almost entirely rattler-free. It should be emphasized that including gravity and body forces obviates the rattler terms in eqs. (1) because then all the grains transmit forces and make cells, with the rattlers making mainly low orders ones. Finally, the theory can be extended to three-dimensions, in which case the main hurdle - cell classification - can be overcome using existing methods Kraynik ; Jordan .

Acknowledgements: C.C.W. acknowledges the Erasmus+ programme of the EU and the Studienstiftung des Deutschen Volkes (German Academic Scholarship Foundation) for funding the visit to the Cavendish Laboratory, where this work was done.