Compliant contact versus rigid contact: A comparison in the context of granular dynamics

We summarize and numerically compare two approaches for modeling and simulating the dynamics of dry granular matter. The first one, the discrete-element method via penalty (DEM-P), is commonly used in the soft matter physics and geomechanics communities; it can be traced back to the work of Cundall and Strack [P. Cundall, Proc. Symp. ISRM, Nancy, France 1, 129 (1971); P. Cundall and O. Strack, Geotechnique 29, 47 (1979)]. The second approach, the discrete-element method via complementarity (DEM-C), considers the grains perfectly rigid and enforces nonpenetration via complementarity conditions; it is commonly used in robotics and computer graphics applications and had two strong promoters in Moreau and Jean [J. J. Moreau, in Nonsmooth Mechanics and Applications, edited by J. J. Moreau and P. D. Panagiotopoulos (Springer, Berlin, 1988), pp. 1–82; J. J. Moreau and M. Jean, Proceedings of the Third Biennial Joint Conference on Engineering Systems and Analysis, Montpellier, France, 1996, pp. 201–208]. The DEM-P and DEM-C are manifestly unlike each other: They use different (i) approaches to model the frictional contact problem, (ii) sets of model parameters to capture the physics of interest, and (iii) classes of numerical methods to solve the differential equations that govern the dynamics of the granular material. Herein, we report numerical results for five experiments: shock wave propagation, cone penetration, direct shear, triaxial loading, and hopper flow, which we use to compare the DEM-P and DEM-C solutions. This exercise helps us reach two conclusions. First, both the DEM-P and DEM-C are predictive, i.e., they predict well the macroscale emergent behavior by capturing the dynamics at the microscale. Second, there are classes of problems for which one of the methods has an advantage. Unlike the DEM-P, the DEM-C cannot capture shock-wave propagation through granular media. However, the DEM-C is proficient at handling arbitrary grain geometries and solves, at large integration step sizes, smaller problems, i.e., containing thousands of elements, very effectively. The DEM-P vs DEM-C comparison is carried out using a public-domain, open-source software package; the models used are available online.

Figure 1

Snapshot of the DEM-P wave propagation in granular media.

Figure 2

The DEM-P wave propagation simulation: balance of total energy $E$ over time. The Coulomb friction coefficient $\mu =0.18$. The total energy was computed as $E=K-U-E_I-W_f-E_A$. Impact, contact, and adhesion are represented using the Hertz, Mindlin, and Derjaguin-Muller-Toporov (DMT) models, respectively [6, 31].

Figure 3

Schematic of a cone drop experiment, with drop height $H$, cone height $L$, and cone width $W$.

Figure 4

Empirical and numerical cone penetration setups showing (a) a cone placed over the settled specimen, (b) a cone penetrating the granular material, (c) fall cones with LVDT connectors, and (d) a view of the assembled apparatus: 1, fall cone and adapter; 2, LVDT; and 3, adjustable vertical stand.

Figure 5

Evolution of the system's kinetic energy (on a lin-$log$ scale) during the first, i.e., settling, stage of the simulation.

Figure 6

Cone depth vs time plots obtained from simulation and laboratory experiments: (a) cone ${30}^{\circ }$ in apex angle, container with a 6-in.-wide diameter, and dense packing; (b) cone ${30}^{\circ }$ in apex angle, container with a 6-in.-wide diameter, and loose packing; (c) cone ${30}^{\circ }$ in apex angle, container with a 4-in.-wide diameter, and loose packing; and (d) cone ${60}^{\circ }$ in apex angle, container with a 6-in.-wide diameter, and dense packing. Black lines indicate experimental data, blue is used for DEM-P results and orange is for DEM-C results. In several cases the blue line is right behind the orange line.

Figure 7

Shear apparatus: 1, sliding portion of the shear box; 2, fixed portion of the shear box; and 3, load cell.

Figure 8

Particle arrangement comparison at the end of the shear test for experiment (top) and simulation (bottom) with granular material (a) randomly packed loosely with an incline angle of ${18}^{\circ }$ and a shearing speed of 1.0 mm/min and (b) hexagonally packed densely with an incline angle of ${24}^{\circ }$ and a shearing speed of 0.5 mm/min. Six particles, marked in red, were selected for analyzes via PIV.

Figure 9

Trajectories for selected particles obtained from simulation and laboratory experiments: (a) particles packed densely and the shear box inclined at ${18}^{\circ }$, (b) particles packed loosely and the shear box inclined at ${18}^{\circ }$, (c) particles packed densely and the shear box inclined at ${24}^{\circ }$, (d) particles packed loosely and the shear box inclined at ${24}^{\circ }$, (e) particles packed densely and the shear box inclined at ${30}^{\circ }$, and (f) particles packed loosely and the shear box inclined at ${30}^{\circ }$. Two plots are provided in (a)–(f); the left and right plots correspond to $V_{\text{sh}}=0.5\text{mm}/\text{min}$ and $V_{\text{sh}}=1.0\text{mm}/\text{min}$, respectively. Two sets of trajectories obtained in simulations are shown. The first set presents the results with shearing velocity 500 times larger than the one used in laboratory ($V_{\text{sh},\text{sim}}=V_{\text{sh}}\times 500$); these simulations are shown with dotted lines, which are always orange for the DEM-C and blue for the DEM-P. The second set of results, shown with solid lines, is obtained for shearing speeds closer to experiment; i.e., the speed-up factors $\alpha$ for $V_{\text{sh}}=0.5$ mm/min and $\beta$ for $V_{\text{sh}}=1.0$ mm/min are both significantly less than 500. Since the DEM-P was slower, the slowest shearing velocities we simulated with were 4 mm/min, which gives $V_{\text{sh},\text{sim}}=V_{sh0.5}\times 8$ and $V_{sh1.0}\times 4$, i.e., $(\alpha =8,\beta =4)$. For the DEM-C, no shearing speed-up was needed in the simulation, i.e., $\alpha =\beta =1$.

Figure 10

Evolution of the system's kinetic energy during the standard triaxial test settling stage.

Figure 11

Granular material simulation snapshot at the end of (a) the settling stage and (b) the standard triaxial test stage.

Figure 12

Laboratory vs simulation results for the standard triaxial test with (a) a monodisperse specimen and (b) a polydisperse specimen.

Figure 13

Flow of prolate ellipsoid through a hopper with $\alpha =3.4$.

Figure 14

Heterogeneous granular flow through a hopper. The material is composed of equal fraction of spheres, ellipsoids ($\alpha =2$), boxes, and cylinders.