Energy minimization and preconditioning in the simulation of athermal granular materials in two dimensions

1.   Introduction

Granular materials are conglomerates of discrete, macroscopic, solid particles, such as sand, soils, pills, etc. They are ubiquitous in many industrial and natural processes. Analytical study is usually difficult due to the complicated nonlinear heterogeneous multi-body interactions in the dense granular system. Numerical methods have remarkable significance in the study of granular materials since most perturbative analytical methods in condense matter physics requires either the low density in the ideal gas limit or a prefect lattice for crystals. On the meanwhile, it is difficult for experimental methods to deal with large number of particles, measure physical quantities in 3D and investigate phase transitions such as jamming phenomena quantitatively [10,25]. Numerical simulations allow access to all detailed information of the particle system, and usually it is straightforward to calculate all the relevant quantities (e.g. fabric, stress, energy) of the system.

However, it remains challenging to develop efficient and robust numerical methods for granular system, especially for large particle systems. Molecular dynamics and energy minimization are two main classes of simulation methods. Discrete element method (soft-particle method) is one of the typical molecular dynamics methods in granular applications [12], where a set of Langevin equations of motion are formulated, with relaxation and random fluctuation terms to quantify the physical relaxation time scale and thermal effect, respectively. Related methods include event-driven method, non-smooth contact dynamics, etc, see references in [11]. On the other hand, energy minimization methods can be applied to the situation of quasi-static deformation, when the applied load changes slowly over time with respect to the inertial forces. One thus can significantly speed up the simulation if time scales are not needed to be fully resolved. Also, the high accuracy of the energy minimization methods provides advantage for the quantitative study of the jamming transition and the scaling law behavior close to the jamming point. Although molecular dynamics method have been extensively used in the simulation of granular materials, the application of energy minimization methods is relatively recent. For example, Luding et. al. have applied trust region methods to study the relaxation and shear of granular system [11]. In this paper, we will study the performance of several typical minimization methods for the granular simulation, such as conjugate gradient methods and quasi-Newton type methods.

Preconditioning is the main bottleneck for the development of efficient and robust numerical algorithms for large scale molecular simulation [1,19], boh in the case of molecular dynamics and energy minimization. General purpose preconditioners are not specifically targeted at large-scale atomistic systems, and are not particularly effective. In this paper, we will also stress on the development of preconditioning techniques for the efficient energy minimization of granular systems.

Remark 1.1. The molecular details of molecular dynamics methods such as soft-particle methods and energy minimization methods are usually not the same. However, it is important to note that the equilibrium physical quantities such as fabric, stress and energy remain similar [11]. Also, there is no temperature in the granular model simulated by energy minimization methods, namely, we are only interested in athermal granular materials.

1.1.   Outline

The paper is organized as follows. In Section 2, we present the physical model for the granular system, as well as the critical jamming transition which is important for the numerical investigation of dense granular system. We present the optimization and preconditioning techniques used for granular simulation in Section 3. The numerical methods are validated and benchmarked in the numerical experiments in Section 4. We conclude in Section 5 and point out some future improvements.

2.   Physical setup

In this section, we will start with a brief introduction for a contact model for granular materials in Sec 2.1 and 2.2, then give the definitions for some macroscopic variables such as stress, elastic modulus and average contact number in Sec 2.3, finally in Sec 2.4 we present the jamming formation which is an important physical phenomena for granular systems. For simplicity, we only consider the two dimensional granular system.

2.1.   Granular configuration

An admissible configuration ${\mathcal{C}}$ of particles is given by a list of center $\vec{x}_i$ and radius $r_i$ of particles. The volume (or packing) fraction $\phi$ is a key parameter in the simulation of granular materials,

where $V$ is the volume of the containing domain (container).

The geometric structure of the granular materials can be characterized by the following fabric tensor,

where $C$ is the collection of all contacts $c$ and

2.2.   Contact model

It is important to note that granular particles can not interact unless they overlap. In other words, there is only repulsive force therefore the system is qualitatively different from the repulsive-attractive forces such as Lennard-Jones [5]. Also for simplicity, we only consider contact models without tangential forces.

Here we will introduce a contact potential of the following form,

where $k_{ij}$ denotes the particle stiffness, $r_i$ and $r_j$ are the radii of particle $i$ and $j$, respectively, $\vec{x}_i$ and $\vec{x}_j$ are the positions of their centers, respectively.

The total potential energy of the system is the sum of interaction energy of all particles.

See Figure 1 for an illustration.

Since the particles interact with each other only if they contact, we can rewrite the potential energy for a single contact $c$ as,

where $\delta_c$ is the overlap associated with the contact $c$, and by equation (4),

In this formula, particle $i$ and $j$ are two particles associated with contact point $c$ (see Figure 1). Therefore, if $\delta_c\geq 0$, the first derivatives, or forces of $E_c$ are,

where $\vec{f}_c$ is the contact force acting on particle $i$ by particle $j$. Furthermore, the second derivatives are,

where $\otimes$ is the tensor product of two vectors and $I$ is the $2\times 2$ identity matrix.

Remark 1. The exponent $\alpha\geq 2$ identifies the nonlinearity of the pair potentials. When $\alpha = 2$, the interaction between particles follows the Hooke's law, which models a repulsive harmonic spring. The force between two overlapped particles linearly depends on the overlap. Hooke's model is one of the simplest contact model, which is still of great importance. From equation (8) and (9), it is clear that the first derivatives of the Hooke potential are continuous at $\delta_c = 0$, while the second derivatives are discontinuous there. When $\alpha > 2$, we obtain the so-called Modlin-Hertzian model or Hertzian model, which is based on the Hertz contact theory [9]. The force between two overlapped particles has a power law nonlinearity. Both the first and second derivatives of the Hertzian potential are continuous if $\alpha >2$. We take $\alpha = 5/2$ in the numerical examples of this paper.

2.3.   Macroscopic variables

The average stress tensor can be calculated by the following formula [2],

where $V$ is the volume of the containing domain and $C$ is the collection of all the contact points in the domain. According to equation (10), the pressure of the system can be defined as[18],

where $d$ is the dimension of the domain. The bulk modulus can then be defined as,

For a simple shear [17] such that the deformation matrix is given by

the shear stress is $-\sigma_{12}$, and the shear modulus can be calculated by:

In addition, the average contact number or coordination number, denoted by $\bar{z}$, is a key geometric parameter in granular simulation. For a granular system consisting of $N$ particles, we can define,

where $|C|$ is the number of all contact points.

We note that there exist particles with zero contacts, namely, the so called "rattlers". In addition, in our 2D frictionless circular disk system, some particles are called "dynamic rattlers" [6] if their numbers of contacts are less than 3, which can lead to mechanical instability. Thus, these rattlers are excluded from the calculation of the coordination number. We denote the number of particles with at least 3 contacts by $N_3$, and the set of contact points for particles with at least 3 contacts by $C_3$. The modified definition of coordination number is,

which will be used in the numerical experiments of this paper.

2.4.   Jamming formation

The jamming transition is an important physical process in granular system [13,18,24], which refers to the critical slowing down of the system due to "overcrowded" particles when the packing fraction $\phi$, defined in (1), increases. Close to the transition point, the geometric constraints prevent the system to explore the phase space, and the system changes its behavior from gas or liquid like to solid like. During the transition, the system has to adjust its configurations more and more in a collective manner instead of via local movements, in order to achieve energy equilibrium. Moreover, This corresponds to the essential difficulty for the numerical simulation which will be discussed later.

In this paper, we adopt the approach by O'Hern et al. [18] to generate a jamming configuration: An initial configuration is sampled with sufficiently low packing fraction ${\phi_{\mathrm{init}}}$ such that there is no overlapping particles after relaxation (energy minimization). The total energy and pressure for the relaxed system with volume fraction ${\phi_{\mathrm{init}}}$ is $0$. Starting from ${\phi_{\mathrm{init}}}$, we repeatedly increase the volume fraction $\phi$ with a small increment ${\delta \phi}$ by enlarging the size of particles accordingly, and then carry out energy relaxation. At some critical fraction ${\phi_{\mathrm{cr}}}$, there will be no free space for particles to explore and unavoidably, they will come into contact. As the system is further compressed, some particles overlap, and the total energy and pressure become nonzero and starts to grow. The system will also posses nonzero bulk modulus (defined in (12)), as the pressure increases upon compression. ${\phi_{\mathrm{cr}}}$ is called the critical volume fraction or jamming point. In 2D bi-disperse systems composed of two kinds of particles with number ratio 1:1 and size ratio 1.4:1, ${\phi_{\mathrm{cr}}}$ is approximately 0.842 [11].

We can adopt a bisection algorithm to increase the accuracy of ${\phi_{\mathrm{cr}}}$. Let ${\phi_{\mathrm{cr}}}^{\mathrm{app}}$ be an approximate critical fraction associated with a configuration with nonzero energy and pressure, the previous unjammed configuration with zero energy and pressure has the volume fraction ${\phi_{\mathrm{cr}}}^{{\mathrm{pre}}}: = {\phi_{\mathrm{cr}}}^{\mathrm{app}} - {\delta \phi}$. We can then compress the unjammed configuration at ${\phi_{\mathrm{cr}}}^{{\mathrm{pre}}}$ by increasing its volume fraction to ${\phi_{\mathrm{m}}} = {\phi_{\mathrm{cr}}}^{{\mathrm{pre}}} + {\delta \phi}/2$, and see if the jamming transition already happens at ${\phi_{\mathrm{m}}}$. This procedure can be repeated until for example, ${\delta \phi}< \varepsilon$, where $\varepsilon$ represents the desired accuracy. We note that this procedure can achieve very high accuracy for the determination of the jamming point, which is not feasible for experimental methods.

The overall procedure is summarized in the Algorithm 1 for the generation of the jamming configuration. Denote by $\min E(\phi, {{\mathcal{C}}_{\mathrm{ref}}})$ the energy minimization in the configuration space ${\mathcal{C}}$ with the packing fraction $\phi$ and the initial configuration ${{\mathcal{C}}_{\mathrm{ref}}}$, which can be solved by certain optimization algorithm in Section 3. The minimizer is denoted by ${\mathcal{C}}_{\mathrm eq} = {\rm arg\;min} E(\phi, {{\mathcal{C}}_{\mathrm{ref}}})$. Here we use a square shaped container $\Omega$ and the periodic boundary condition for the granular system.

3.   Numerical method

The minimization of the total potential energy (5) relies on the state-of-the-art optimization and precondition techniques which will be presented in this section. The performance of those methods is contingent upon the particular applications. In particular, the construction of a good preconditioner, can be regarded as "an art rather than a science" [23].

3.1.   Energy minimization

In this paper, we use the limited-memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) method [15] and Fletcher-Reeves conjugate gradient (FR-CG) method [16] for the energy minimization of granular system. We denote the total energy of a system as $f(x)$ where $x$ is the positions of all particles. If $x^{({k})}$ is the iterate at step $k$, we denote the energy at step $k$ by $f^{({k})} = f(x^{({k})})$, and the gradient at step $k$ by $g^{({k})} = \nabla f^{({k})}$.

3.1.1.   L-BFGS method

L-BFGS is one type of the quasi-Newton methods, namely, it utilizes an approximation of the inverse Hessian matrix to generate the search direction. Instead of the full dense $n\times n$ approximations to the inverse Hessian ($n$ is the number of variables, which is $Nd$ with $N$ number of particles and $d$ the dimension in our problem), L-BFGS uses a low rank approximation with only a few vectors of length $n$ to represent the approximation implicitly.

We denote the approximate inverse Hessian at step $k$ by $B^{(k)}$. Assuming we have stored the last $m$ updates of the form $s^{({k})} = x^{({k+1})}-x^{({k})}$ and $y^{({k})} = g^{({k+1})}-g^{({k})}$, we define $\rho^{({k})} = \frac{1}{y^{(k),T}s^{({k})}}$. The L-BFGS search direction $p^{{({k})}} ( = B^{(k)}\nabla f^{({k})})$ can be obtained using the following algorithm.

The matrix $B^{(k)}_0$ in the boxed step in the algorithm is an rough approximation to the inverse of the exact Hessian, and an effective choice is:

This choice ensures that the search direction is well scaled and therefore the unit step length is accepted in most iterations. Besides, the diagonal matrix makes it much simpler to compute the multiplication $r = B^{(k)}_0q$.

3.1.2.   FR-CG method

We use the Fletcher-Reeves variant [7] of the nonlinear conjugate gradient (CG) method [20]. As an extension of the linear conjugate gradient method [8], the search direction $p^{(k)}$ at $k$th step of FR-CG method is defined by,

3.1.3.   Line search

Once we calculate the new search direction using either L-BFGS in Algorithm 2 or FR-CG in Algorithm 3, the next iteration is given by

where $\alpha^{({k})}$ is chosen by a line search method which satisfies Wolfe conditions or strong Wolfe conditions [16], such that the update is stable.

3.2.   Preconditioning

Preconditioning is important for the efficiency of the minimization algorithms, especially for large scale problems. We have to choose the preconidtioner matrix $P^{({k})}$ such that it is similar to the exact Hessian, and at the same time easy to invert.

3.2.1.   Preconditioned minimization algorithm

For L-BFGS algorithm 2, we replace the boxed step in the algorithm with:

in order to obtain a preconditioned L-BFGS algorithm.

For FR-CG algorithm 3, we have the following preconditioned version for $p^{{({k})}}$.

The equation (16) or Step 1 in Algorithm 4 is equivalent to solve a linear system:

An alternative point of view for the preconditioning is to make a change of coordinates: ${\widetilde{x}} = P^{(k),1/2}x$. Considering the function $F({\widetilde{x}}) = f(P^{(k),-1/2}{\widetilde{x}})$, we have

Applying the L-BFGS algorithm 2 or FR-CG algorithm 3 to optimize the transformed function $F(\tilde{x})$, we can obtain the preconditioned version of the corresponding algorithm.

3.2.2.   Construction of preconditioners for large scale molecular system

Preconditioning is well established in numerical linear algebra and numerical PDE problems [23]. However, in many applications general purpose preconditioners do not work particularly well, and specifically designed preconditioned are preferable. For large scale atomic/molecular simulations, Packwood et. al. proposed the so-called "universal preconditioner" $P$ [19]. They have successfully applied this preconditioner to atomic simulation and electronic structure calculation of crystalline system of the size from hundreds to $10^4$ atoms, with a speedup of order 1 or 2 [14]. In this paper, we will test the performance of those universal preconditioners for granular systems.

The universal preconditioner $P$ is defined via the quadratic form

Here $r_{{\mathrm cut}}$ is the cutoff radius, $\mu$ is the energy scale, $A$ is a parameter which should be large enough to ensure that the nearest neighbor interactions dominant. The numerical experiments in [19] and in our numerical results indicate that as long as $A$ is of order 1, it does not significantly affect the performance of the preconditioner. In particular, if we take $A = 0$, we choose the preconditioner matrix $P$ as the stabilized adjacency matrix of the particles [19], namely,

where $d_{ij}$ is the distance between particle $i$ and $j$, $C_{{\mathrm stab}}$ is a stabilization term to make sure the matrix is strictly positive definite, we choose $C_{{\mathrm stab}} = 0.1$ in our application.

In our model problem, the bi-disperse granular system is composed of two types of particles, whose radii are denoted by $r_{\min}$ and $r_{\max}$. We take the cutoff radius $r_{{\mathrm cut}} > 2r_{\max}$ so that all the possibly overlapping particles are covered. For example, we can take $r_{{\mathrm cut}} = 2.2r_{\max}$. The preconditioner matrix will be recalculated when the maximum displacement of particles is beyond $r_{cut}/4$, therefore it will not be recalculated frequently. The energy scale $\mu$ is chosen to make sure that the precondition matrix is of the same order as the real hessian, so we may choose the unit step-length when the L-BFGS method is applied. One way to obtain $\mu$ is,

where $P_{\mu = 1}$ is the matrix when $\mu = 1$, and $\nu$ is a test perturbation with the following form

where $L$ is the lengths of the periodic cell and $M$ is a constant, here we choose $M = 10^{-3}r_{cut}$.

Since we need to solve linear system (17) with respect to $P$, when the particle number is small, we use direct method and rearrange the index of grains by sparse reverse Cuthill-McKee ordering method [4], the band width of $P$ can be reduced so that the equation (17) can be easily solved. If the number of particles is large, iterative solvers such as preconditioned CG or AMG [3] can be utilized.

4.   Numerical results

In our numerical experiments, we use the bi-disperse granular system which consists of two types of particles with number ratio 1:1 and size ratio $r_{\max}/r_{\min} = 1.4$, so the system will not crystallize [21]. The particles are contained in a unit square domain $\Omega = [0,1]^2$. We use the periodic boundary condition in order to reduce the boundary layer effect.

In the numerical experiments, we test the performance of preconditioned L-BFGS method and preconditioned FR-CG method, by running several benchmark simulations: the jamming formation, the scaling law behavior for the jammed configuration, and the shear deformation of the over-jammed states. In those different scenarios, we will show that the proposed preconditioned algorithm has better performance compared to their un-preconditioned counterparts, and preconditioned L-BFGS is currently the method of choice.

4.1.   Jamming formation

We first run Algorithm 1 to generate the jamming configuration, starting from a relaxed configuration with volume fraction $\phi< {\phi_{\mathrm{cr}}}$, where ${\phi_{\mathrm{cr}}}$ is the critical volume fraction (or jamming fraction). We let the particles grow with small increment $\delta\phi$, then we use either the preconditioned or un-preconditioned version of Algorithm 2 and Algorithms 3, 4 as the energy minimizer in Algorithm 1.

In Figure 2, the residual norm of each iteration during a relaxation process is plotted. We apply four methods (L-BFGS, preconditioned L-BFGS, FR-CG, Preconditioned FR-CG) to minimize the energy of an unjamming configuration. A high accuracy tolerance is adopted here when an accurate jamming point is pursued. The figure shows that: L-BFGS is better than FR-CG, preconditioned L-BFGS method has the best convergence curve compared to other methods, and preconditioner provides a factor of two speed up for L-BFGS.

We also notice that convergence curves for all methods are extremely rough and have long asymptotic regimes, that indicates the granular system may have a complex and rough energy landscape.

Figure 3 shows the evolution of the granular system up to the jamming point. The volume fraction of the initial configuration is small, therefore particles can move around in order to achieve a zero energy configuration. When the volume fraction increases, the particles start to contact with each other, but initially it is still possible to achieve a zero energy configuration through energy relaxation (minimization). As the free space becomes less and less, at a critical volume fraction, the particles are forced to contact and overlap, and the potential energy of the system becomes non-zero.

4.2.   Scaling law

Once the jamming configuration ${{\mathcal{C}}_{\mathrm{jam}}}$ is obtained, we continue increase the volume fraction quasi-statically. It has been proposed by physicists that granular systems have the so-called scaling law behavior in the vicinity of the jamming point [18]. Namely, macroscopic physical quantities such as pressure $p$, shear modulus $G$ and average contact number $z$ defined in equation (14), vary as a power of volume fraction $\phi-{\phi_{\mathrm{cr}}}$. We refer to [18] for more detailed discussions of the physical interpretations of scaling law behavior near critical points. In this paper, we numerically justify the scaling laws and calculate the critical exponents using our energy minimization techniques. Interestingly, we observe that the iteration numbers of the energy minimization methods also has similar scaling law behavior.

In Figure 4, we plot the log-log curve of pressure $p$ vs $\phi-{\phi_{\mathrm{cr}}}$. It is clear that $p$ follows a power law with respect to $\phi-{\phi_{\mathrm{cr}}}$ close to the jamming point, and the critical exponent depends on the nature of interactions. For harmonic $(\alpha = 2)$ and Hertzian $(\alpha = 2.5)$ potentials, the critical exponent are approximately $1.0$ and $1.5$, respectively. In addition, this scaling law between $p$ and $\phi-{\phi_{\mathrm{cr}}}$ is robust for different initial configurations (they may have different ${\phi_{\mathrm{cr}}}$).

If we exert small simple shear strain $\gamma$ to these configurations, we can compute their shear modulus using equation (13). In Figure 5, we plot shear modulus for different volume fractions. It shows that the shear modulus $G$ of the granular systems also follows a similar power scaling law with respect to $\phi-{\phi_{\mathrm{cr}}}$ but with different exponents, which is approximately $0.5$ when $\alpha = 2.0$, and $1.0$ when $\alpha = 2.5$.

The modified average contact number $z$ defined in (14) also exhibits a scaling law in regards to volume fraction difference $\phi-{\phi_{\mathrm{cr}}}$. In Figure 6, the lines in the log-log plots have similar slopes close to $0.5$ for both $\alpha = 2.0$ and $\alpha = 2.5$. The critical average contact number $z_0 = 4$ is the average contact number when a system reaches exactly the jamming state. The average contact number of a system takes a leap from zero to 4 when jamming transition occurs.

In Figure 7, we plot the iteration number vs. volume fraction difference $\phi - {\phi_{\mathrm{cr}}}$ for the preconditioned L-BFGS and L-BFGS methods. Starting from certain jamming configuration, we compress the granular system by increasing the volume fraction $\phi$. Then for each relaxed over-jammed configuration, we exert small shear strains from $10^{-6}$ to $10^{-5}$ quasi-statically, in 10 steps. The iteration numbers exhibit "random" behavior, especially close to the jamming point. We plot both the average and the variance (as error bars) of the iteration number. We observe that the iteration number of each numerical method also follows a scaling law, which ranges from $10^3-10^4$ close to the jamming point ($\phi-{\phi_{\mathrm{cr}}} = 10^{-6}$), to about hundreds when $\phi-{\phi_{\mathrm{cr}}} = 10^{-2}$. The result shows that the preconditioned algorithm performs better, and also it is more robust to rounding error as the variance of the iteration number for the preconditioned algorithm is also smaller away from the jamming point.

In Table 1, we illustrate the performance of four optimization methods: L-BFGS, preconditioned L-BFGS, FR-CG, and preconditioned FR-CG in three different situations given an over-jammed configuration with 4096 particles and volume fraction $\phi = 0.8438$.

● Case 1, fixing the positions of the particles, and enlarging them slightly with ratio $1+\delta, \delta = 5\times 10^{-5}$;

● Case 2, applying a small shear strain ($\gamma = 10^{-4}$) to the configuration;

● Case 3, exerting a small ($\mathcal{O}(10^{-4})$) random perturbation to the configuration.

The comparison of iteration number and computational time in Table 1 clearly demonstrates that the preconditioner can improve the efficiency of the energy minimization methods. Also, it seems we can simply take the parameter $A = 0$ (19) in the preconditioner.

4.3.   Shear test

In the shear test, we deform the over-jammed configuration with pure shear (stretching in x direction and compressing in y direction, while keeping the area unchanged). The deformation matrix is given by,

The stress vs. strain and fabric vs. strain curve for the pure shear simulation is shown in Figure 8. Notice that we use increments of different sizes, $4\times 10^{-4}$, $2\times 10^{-4}$, and $10^{-4}$ to reach a final 2% shear strain, we still obtain qualitatively (quantitatively for small strain) similar curves.

5.   Conclusion

In this paper, we introduce the energy minimization techniques for the efficient simulation of dense athermal granular systems in two dimensions. Preconditioners are used to enhance the performance of the simulation. We carry out numerical experiments for some typical scenarios of granular simulation in order to validate the numerical methods, which include the jamming formation, scaling law phenomena close to the jamming point, and deformation of over-jammed states. Speedup of the preconditioned minimization methods are observed.

This work opens avenue for several possible improvements: First of all, we are going to extend the current work to three dimensions, which is physically more relevant to real applications, and still difficult for experimental studies [10,25].

Secondly, from the practical point of view, we will optimize the implementation of preconditioners using, for example, AMG [3] and other types of state-of-the-art techniques, especially for the three dimensional case.

Last but not least, we observe that the iteration number of our current preconditioned algorithms follows a scaling law behavior close to the jamming point, which is a fundamental bottleneck for the energy minimization approach. It remains open whether one can find a preconditioner which is robust close to the jamming point. Usually, an efficient preconditioner takes account of the long wavelength modes of the system. However, close to the phase transition point, more and more high frequency information might be needed. The development of numerical techniques in this direction relies on the understanding of the physical origin of the jamming transition.

Acknowledgments

We thank Jie Zhang and Zhaohui Jin from Shanghai Jiao Tong University for stimulating discussions on granular modeling, experiments and simulation. We thank Christoph Ortner from University of Warwick and Mingjie Liao from Shanghai Jiao Tong University for stimulating discussions on numerical methods of atomic simulation and their efficient implementations.