DEM simulations of polydisperse media: efficient contact detection applied to investigate the quasi-static limit

Discrete element modeling (DEM) of polydisperse granular materials is significantly more computationally expensive than modeling of monodisperse materials as a larger number of particles are required to obtain a representative elementary volume, and standard contact detection algorithms become progressively less efficient with polydispersity. This paper presents modified contact detection and inter-processor communication schemes implemented in LAMMPS which account for particles of different sizes separately, greatly improving efficiency. This new scheme is applied to the inertial number (I), which quantifies the ratio of inertial to confining forces. This has been used to identify the quasi-static limit for shearing of granular materials, which is often taken to be I=10-3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ I = 10^{ - 3} $$\end{document}. However, the expression for the inertial number contains a particle diameter term and therefore it is unclear how to apply this for polydisperse media. Results of DEM shearing tests on polydisperse granular media are presented in order to determine whether I\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ I $$\end{document} provides a unique quasi-static limit regardless of polydispersity and which particle diameter term should be used to calculate I\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ I $$\end{document}. The results show that the commonly used value of I=10-3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ I = 10^{ - 3} $$\end{document} can successfully locate the quasi-static limit for monodisperse media but not for polydisperse media, for which significant variations of macroscopic stress ratio and microscopic force and contact networks are apparent down to at least I=10-6\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ I = 10^{ - 6} $$\end{document}. The quasi-static limit could not be conclusively determined for the polydisperse samples. Based on these results, the quasi-staticity of polydisperse samples should not be inferred from a low inertial number as currently formulated, irrespective of the particle diameter used in its calculation.


Introduction
Polydisperse granular materials (i.e., those containing a range of particle sizes) occur in many physical and industrial settings, such as geomaterials, avalanches and landslides, crushing of mining ores and food processing. Often the range of particle sizes in such systems covers several orders of magnitude. For example, the granular filter material used to construct the Bennett Dam in Canada contains particles ranging from 0.08 to 75 mm [1]. Such a range of length scales presents significant computation challenges to the discrete element method (DEM), typically used to model such systems. These challenges reflect the fact that in polydisperse systems (i) a larger number of particles are required to obtain a representative elementary volume (REV) than for a monodisperse system and (ii) the standard contact detection algorithms used in such modeling can become progressively less effective with increasing particle size ratio.
While growing computational power has allowed effective investigation of increasingly polydisperse systems [2][3][4], and algorithmic enhancements in contact detection [5,6] have improved the efficiency of such simulations, there remain challenges. This remains particularly true for simulations requiring long timescales. Many physical processes and standard laboratory tests such as geomechanical element testing involve extremely low strain rates imposed over a long time period, and it can generally be assumed that quasi-static shearing occurs. These conditions would be impractical to replicate in a DEM simulation of reasonable computational cost, so the simulated strain rates are artificially increased by orders of magnitude. Correspondence between the simulations and reality is maintained by loading the granular material quasi-statically, i.e., the loading occurs sufficiently slowly that inertial effects can be neglected. In order to identify the boundary between the quasi-static and inertial or dense-flow regimes, the dimensionless inertial number I εd ρ p is used, whereε is the shear rate, d is particle diameter, ρ is particle density and p is the mean confining stress [2,3]. The inertial number represents the ratio of inertial to confining forces, and as I → 0, the flow regime tends to the quasi-static limit. Radjai [4] states that "For a confining pressure p (counted positive for compressive stresses) and particles of average diameter d, the contact forces of static origin are of the order of f s pd 2 … At the same time, for a shear strain rate˙ q , the timescale of the flow is and thus the order of magnitude of the impulsive forces is given by the momentum per unit time f i md˙ q / t, where m is the average particle mass. In the quasi-static limit, the condition f s f i implies I ≡˙ q m pd 1." Andreotti et al. [5] explain that I can be interpreted as the ratio between two timescales: I representing the rate of microscopic rearrangements of particles subject to a pressure p, and t macro 1 ε , representing the macroscopic shear rate. In the quasi-static regime, macroscopic rearrangements can be considered to occur very slowly in comparison with the microscopic rearrangements.
Based on the empirical assessment of 2D DEM simulation results of plane shear tests, da Cruz et al. [2] set the practical limit of the quasi-static regime at I ≤ 10 −3 . Considering conditions at the critical state in 3D DEM simulations of geotechnical element testing, Perez et al. [3,6] found the limit at I ≤ 7.9 × 10 −5 .
The above work has been influential in improving the quality of DEM simulations for granular materials under quasi-static conditions in that it is now relatively common to set shear rates so that I ≤ 10 −3 or I ≤ 10 −4 . However, the use of a single particle diameter, d, in the calculation of I suggests a monodisperse material. In reality, many granular materials, including most geomaterials, are polydisperse. An ideal definition of inertial number would be able to identify a unique quasi-static limit regardless of the particle size distribution of the material under shear. The selection of an appropriate diameter term for use in the calculation of inertial number is required to define such an inertial number. Apart from the work of Rognon et al. [7], who proposed a packing fraction-weighted inertial number to account for granular flows involving disks of a different diameter to account for segregation, there has been very little work to examine the effectiveness of inertial number in defining the quasi-static limit for polydisperse materials.
To provide access to the low strain rates and long timescales required to address the question of where the effective inertial regime lies in polydisperse systems, this paper first addresses an improved contact detection method implemented in the popular molecular dynamics code LAMMPS [8]. A series of DEM triaxial compression tests is then carried out at varying shear rates on materials with varying degrees of polydispersity using the new contact detection method. Analysis of a selection of preliminary results allows the validity of the previously proposed limits for quasi-static behavior to be examined using various diameter terms to calculate the inertial number.

Modified DEM model
DEM simulations were carried out using a modified version of the open-source DEM code Granular LAMMPS [8]. In order to improve efficiency when simulating highly polydisperse materials, the contact detection and inter-processor communication algorithms in LAMMPS were modified from the existing link-cell method [8] to a new method termed the hierarchical stencil method, which is conceptually similar to the hierarchical grid Method used in MercuryDPM [9,10] but utilizes the existing LAMMPS stencil capabilities developed by in't Veld et al. [11]. A brief outline of the modification is given here; full details of the implementation and parametric studies are available at [12].
The existing LAMMPS contact detection is a combination of the widely used Verlet neighbor list and link-cell methods [8]. In DEM simulations, Verlet neighbor lists store all pairs of particles which are within a distance 2r skin of each other, where r skin is the "skin distance," as defined in Fig. 1. The additional skin distance means the neighbor list must be constructed intermittently (e.g., when any particle has moved a distance of r skin /2 since the last rebuild). At each intermediate timestep, only the particle pairs on the neighbor list are checked for contact; where contact exists, the force is calculated. To avoid brute-force construction of the neighbor list, the link-cell method is used. A regular grid of cells is overlaid on the DEM domain, and, for each particle, a subset of link cells are searched to create the neighbor list. For example, in Fig. 1 the link-cell length is r cell r max + r skin . Considering the green particle in Fig. 1, the neighbor list is constructed by checking the "home" link cell plus surrounding link cells within 2 r c . The list of cells to be checked is stored in a pre-computed stencil [11].
LAMMPS implements a standard approach of domain decomposition with message passing via the message passing interface (MPI) in parallel. This means that particles are "owned" by a given MPI task dependent on their posi-

Fig. 1 Schematic showing link-cell contact detection in LAMMPS,
where r max is the maximum particle radius, and r c is the cell dimension. r c r max + r skin where r skin is the skin distance Particle information within the halo of dimension r halo must be communicated to the processor subdomain at every timestep tion. In order that all relevant interactions may be located, a given local domain must obtain information on particles in adjoining regions. Communication of ghost particles within a "halo" of cells with dimension r halo 2r c is performed to fulfill this requirement, as shown in Fig. 2.
The Verlet/link-cell method is highly efficient for monodisperse packings of particles. However, as the pack-ings considered become increasingly polydisperse, the efficiency of the method reduces [12]. Consider a granular system with two types of particle having radii r s and r l (for small and large). This introduces three different interaction cutoff ranges: r ss c , r sl c and r ll c . In principle, one should choose the cell width to be r cell 2r ll c to ensure all large-large interactions are captured. However, the resulting cell size will necessarily drag into the search very many small-large and small-small pairs well beyond their respective cutoffs. This is inefficient and becomes more inefficient as the ratio r l /r s increases. Similarly, the communication halo will be of dimension r halo 2r ll c . Therefore, as the link-cell and halo sizes are both based on the largest particle, for polydisperse packings many more particle pairs must be considered in neighbor list construction and inter-processor communication than for monodisperse packings.
The new hierarchical stencil method overcomes this limitation as follows: • Particle types are allocated based on particle radius; • Cell lists are instantiated for each particle type. The sizes of the link cells are based on the largest particle of each type in a similar way to Ogarko and Luding [9]. For example, for a bidisperse system, two particle types and two cell lists with sizes r s c and r l c are instantiated. • Interactions between particles of the same type are identified using a stencil within the appropriate cell list as shown schematically in Fig. 3 for a bidisperse system. • For interactions between two particles of different types, particle i is located within the cell list of the j-type particles. An appropriate stencil is then used to perform the neighbor list construction in the j-type cell list. The most efficient way to locate particles is using a one-way search which identifies potential small-large pairs by considering only small particles and using the large particle cell list to search for interactions, and using a symmetric stencil in the large cell list (Fig. 3b) to examine large neighbors [12]. • In't Veld et al. [11] improved the existing LAMMPS interprocessor communication by introducing multiple halos for different interaction types. For a bidisperse system, the halo width is r sl halo for small particles and r ll halo for large particles, where, for example, r sl halo r s c + r l c as shown in Fig. 4. This is sufficient to allow identification of all potential pairs. Potential small-small pairs may be located on the basis of r ss halo . In addition, a significant efficiency saving can be made as all potential small-large pairs are located by examining owned small particles, meaning no small ghost particles beyond r ss c are required. The ghost cutoff distance for large particles is unchanged at r ll halo , and potential large-large pairs are identified as before.
• The discussion thus far has focused on bidisperse systems.
Generalization of the scheme to polydisperse systems is straightforward. A number of cell lists of varying size are c large-large interactions selected, and particle i is assigned a type based on the smallest cell for which r i + r skin ≤ r c . The number and size of cell lists must then be selected. Previous studies [13] suggest this is strongly dependent on the particle size distribution for the problem at hand, and no general law is available to decide without testing. However, for the continuous polydisperse systems simulated here it was found that two or three cell lists with a logarithmic size spacing were optimal [12]. This is in contrast to the findings of Krijgsman et al. [13] for the hierarchical grid method, highlighting that although the two schemes are conceptually similar, important differences in implementation exist.
More detail on the classes of the C ++ implementation in LAMMPS is given in [12]. The implementation was vali- Halos of different dimensions are adopted depending on the interaction type (i.e., large-large, small-large or small-small) dated using the analytical solution developed for the failure stress ratios in a face-centered cubic assembly of uniform rigid spheres [14] in which multiple particle types were assigned. A further validation was carried out by comparing the results of triaxial compression simulations using 74,504 particles to the existing link-cell contact detection schemes. The variation in coordination number and stress ratio at the critical state were found to be a fraction of a percent, representing rounding errors accumulated due to the neighbor lists being assessed in a different order with the different schemes.
The speedup of the hierarchical stencil method over the link-cell method improves with increasing size ratio. For size ratios of r max /r min 10, speedups were at least 10, and for r max /r min 100, speedups of up to 400 versus the existing link-cell method without communication improvements were obtained [12]. The hierarchical stencil also scales well to at least 768 processors, with scaling being greatly improved by the inter-processor communication improvements [12].

DEM simulations
A total of 76 DEM simulations of constant mean stress triaxial tests were carried out. Seven different polydisperse particle size distributions (PSDs) were simulated, as shown in Fig. 5. Samples of series "A" have an equal volume of particles per log diameter bin, whereas samples of series "B" have an equal volume of particles per linear diameter bin. Series A therefore have relatively more fine particles. The  Table 1. These were selected on a trial-and-error basis, taking care to ensure that an REV was achieved for each sample so that the sample responses with respect to I are meaningful. Unfortunately, no clear relationship can be established between grading and number of particles for an REV. A conservative timestep of 7.5 × 10 −8 s was used in all simulations, calculated using t 0.1 m min K max where m min is the minimum particle mass and K max is the maximum contact stiffness [15] calculated using a 2% particle overlap (actual overlaps in the simulations at no point exceeded 1%).
In each test, the particles were initially generated in a random, non-touching cloud, before being isotropically compressed to p 100 kPa. A simplified Hertz-Mindlin contact model was used with shear modulus G 29 GPa and Poisson's ratio ν 0.2. The initial interparticle friction coefficient of μ 0.15 was used during compression to create an initially dense packing configuration. Following isotropic compression, the friction was set to μ 0.3 and the sample allowed to equilibrate. During shearing, the mean normal stress was maintained at p 100 kPa to allow a constant value of I to be maintained [3]. For each PSD, a series of tests were carried out in which the axial shear rate,ε 1 , was varied to impose different values of inertial number calculated using the maximum particle diameter, I dmax ε 1 d max ρ p , ranging from I dmax 5 × 10 −3 to I dmax 5 × 10 −6 for samples with χ d max /d min 1.2 to 5 and I dmax 5 × 10 −3 to I dmax 1×10 −6 for samples with χ 10 and 20. I dmax was selected as the default inertial number as it gives the largest value of I. I dmax 1 × 10 −6 was proved to be the slowest simulation which could be practically carried out: To shear to ε 1 2% at this rate, sample A20 required around 480 h using 180 cores and B20 required 288 h using 72 cores on the Cirrus HPC facility (http://www.cirrus.ac.uk/). To reduce I dmax by a further order of magnitude would have been computationally infeasible.

Results and discussion
Plots of axial strain against stress ratio η q p , where q is the deviatoric stress, volumetric strain ε v and mechanical coordination number Z mech (the average number of contacts per stress-transmitting particle) for simulations with I dmax 1 × 10 −3 are shown in Figs. 6, 7 and 8. These strains are sufficient to allow the effect of I on material behavior to be determined in what Roux [16] called Regime 2, in which particle rearrangements control the macroscale quasi-static response. At this fixed inertial number, η and |ε v | increase, while Z mech decreases with increasing χ . Figure 9 shows the stress ratio at ε 1 2% with varying I dmax for all samples. The stress ratios are normalized by their val- ues at I dmax 1×10 −3 , which are shown in Fig. 6. The most uniform sample, A1.2, shows approximately constant values at I dmax ≤ 1 × 10 −3 . For all other samples, the stress ratio reduces with inertial number with no sign of a plateau in values at I dmax 1 × 10 −6 , most significantly for sample B20 for which η 0.846η (1e−3) . Interestingly, the trends do not exactly follow χ ; most notably, A10 shows more variation with I dmax than A20. Figure 10 shows normalized values of solid packing fraction φ at ε 1 2%. Apart from the nearmonodisperse A1.2, all samples show an increase in packing fraction as the inertial number reduces, although the magnitude of this change is less significant than for the stress ratio. In contrast, mechanical coordination number (Fig. 11) shows a similar trend regardless of the particle size distribu- Further insight into changes in the stress-transmitting fabric of the sample showing the greatest variation in η with I dmax , B20, is given in Figs. 12a, b and 13a, which, respectively, show the probability density functions of normal contact force and the relative frequency distribution of connectivity (number of stress-transmitting contacts per particle). In Fig. 12a, b it can be seen that the force network tends to become more inhomogeneous as the inertial number increases (resembling the increasing force inhomogeneity found by Voivret et al. [17] in 2D and Mutabaruka et al. [18] in 3D with increasing polydispersity). This suggests that at higher inertial numbers the deviatoric stress-transmitting strong-force network is more dominant, which explains the higher stress ratios at higher inertial numbers. Despite the relatively small changes in mechanical coordination number (Fig. 11), the number of contacts per particle shows significant differences between samples sheared with different inertial numbers (Fig. 13a). The more slowly a sample is sheared, the fewer relatively unstable particles with C 2 or 3 or highly connected particles with C ≥ 18 are present. How- Fig. 11 Effect of I dmax on mechanical coordination number at ε 1 2%. Results are normalized by the response at I dmax 1 × 10 −3 : a series A; b series B ever, there are more particles with 4 ≤ C ≤ 17 when inertial numbers are low. In contrast, the near-monodisperse sample A1.2 has almost indistinguishable force and contact distributions for all samples with I dmax ≤ 1 × 10 −3 as shown in Figs. 12c and 13b. Considering both macro-and microscale results, it can be concluded that a true quasi-static limit has not been reached for the polydisperse samples with χ ≥ 5. For frictional particles, the minimum number of contacts for static mechanical stability is 4 [19], and therefore, the number of particles with four or more contacts can be taken as a measure of quasi-staticity [3,20]. This was termed the nonrattler fraction f NR by Bi et al. [21] and is plotted against inertial number normalized by values at I dmax 10 −3 in Fig. 14a and as raw values in Fig. 14b. Two features to note are: (i) More polydisperse samples have a much lower f NR (i.e., a greater proportion of rattlers) and (ii) f NR reduces with I dmax for the more polydisperse samples. Shen and Sankaran [22] demonstrated that at higher strain rates the coordination number reduces, but the size of groups of interconnected particles (analogous to the non-rattlers) increases, similar to the trend seen here. Large numbers of rattlers will naturally be present in all highly polydisperse materials. It is possible that at inertial numbers around the previously defined quasi-static (monodisperse) limit (I 10 −3 ) these rattlers are more able to join and stabilize buckling force chains [23, 24] than at either higher or lower inertial numbers. This would account for the higher stress ratio and packing fraction at inertial numbers close to the monodisperse limit. These rattlers would be the smaller particles, which would be "captured" by the larger particles upon force chain buckling [25]. Interestingly, at I dmax 5×10 −3 (above the usual definition for the quasistatic limit) the non-rattler fraction and stress ratio are both higher, suggesting that this rattler "capturing" mechanism is mainly found below the monodisperse quasi-static limit. However, the relationship between rattlers and quasi-staticity requires further study. Figure 15 presents the variation of normalized stress ratio for the series B tests with inertial number where two alternative definitions of inertial number are used: (i) I d50 ε 1 d 50 ρ p , Fig. 13 Relative frequency plot of connectivity, C: a sample B20, contacts C ≤ 6 (note linear y-scale); b sample B20, contacts C ≥ 6 (note log y-scale); c sample A1.2, contacts C ≤ 6 (b) (a) 10

Alternative definitions of inertial number
, N mech is the number of particles with two or more contacts, d i is the diameter of particle i and V p,i is the volume of particle i. Id takes a form similar to that proposed by Rognon et al. [7] for bidisperse granular flows. For both I d50 and Id, there is a similar trend of reducing η with inertial number, but the minimum inertial numbers are lower for I d50 and Id than for I dmax . As the quasi-static limit has not been reached for samples with χ ≥ 5, the most effective inertial number for determining this limit cannot be established. As explained in the Introduction, the fundamental concept of inertial number is a ratio between the impulsive forces and the contact forces of static origin [4]. For a polydisperse granular material, the largest possible inertial number, as currently defined, requires maximizing the order of magnitude of the impulsive forces by using the largest particle diameter and mass in its calculation, while minimizing the contact forces of static origin by using the smallest particle diameter. In that "worst possible" case, the inertial number would be χ I dmax . Therefore, the current definition of inertial number does not permit differences in Fig. 15 Effect of alternative inertial number definitions on stress ratio at ε 1 2%: a I d50 ; b Id inertial number of more than three orders of magnitude compared to the uniform case, irrespective of the definition of particle diameter adopted. Hence, the inertial number, as currently defined, is not appropriate for locating the quasi-static limit for polydisperse granular materials.

Conclusions
Particulate simulations of continuously polydisperse granular materials with χ d max /d min 1.2 to 20 were carried out using a DEM code which was modified for increased efficiency with polydisperse media. This was achieved by introducing a hierarchy of cell lists and improved interprocessor communication for particles of a different diameter. In order to investigate whether the quasi-static limit is the same for granular materials regardless of their particle size distribution, the polydisperse samples were sheared under triaxial compression to ε 1 2% with inertial numbers calculated using the maximum particle diameter ranging from I dmax 5 × 10 −3 to I dmax 1 × 10 −6 . From the results, the following conclusions can be drawn: • For a near-monodisperse particle size distribution (χ 1.2), the quasi-static limit was found at approximately I dmax ≤ 1 × 10 −3 , in agreement with previous studies [2]. • For the more polydisperse distributions, the quasi-static limit was not found even at I dmax 1 × 10 −6 and, in general, more polydisperse distributions showed a greater reduction in stress ratio and more homogeneous force and contact networks with a reduction in inertial number.
• More polydisperse distributions have more rattlers, and the proportion of rattlers increases as inertial number reduces below I dmax 1 × 10 −3 . These rattlers may be less likely to join and stabilize force chains at low inertial numbers, leading to a lower stress ratio. • Definitions of inertial number using alternative diameter definitions, for example the median particle diameter, I d50 , or a volume-weighted diameter, Id , are also unable to determine a unique quasi-static limit regardless of particle size distribution.
As currently defined, the inertial number is not appropriate for locating the quasi-static limit for polydisperse granular materials. Further work is required to determine where the quasi-static limit lies for polydisperse media and to establish whether the inertial number could be somehow adapted to find this limit accounting for polydispersity. As both computational resources increase in power and further algorithmic improvements can be identified, it is hoped that future work will be able to access more highly polydisperse systems at yet smaller inertial numbers. Such simulations should be able to identify more exactly where the limiting inertial number lies.